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Frozen  Orbits — Near  Constant  or  Beneficially  Varying  Orbi'al  Parameters 
Thesis  directed  by  Professor  Robert  D.  Culp 

If  a  satellite  were  experiencing  pure  Keplerian  motion  its  orbital  elements 
would  remain  constant.  In  actuality,  each  orbiting  body  is  acted  upon  by  various 
perturbing  forces.  Consequently,  the  orbital  elements  are  continuously  subject  to 
change.  The  forces  generating  deviation  from  basic  two-body  motion  consist  of,  but 
are  not  limited  to:  (1)  the  Earth’s  oblateness,  (2)  atmospheric  drag,  (3)  lunar-solar 
gravitation,  and  (4)  solar  radiation  pressure. 

Whereas  the  magnitudes  of  such  perturbing  forces  are  small  and  varying 
compared  to  the  gravity  effects  of  the  Earth,  they  are  persistent  and  cause  the  resulting 
orbit  to  depart  from  its  Keplerian  counterpart.  The  variation  of  the  five  orbital  elements 
which  delineate  the  size,  shape,  and  orientation  of  an  elliptic  orbit,  can  be  determined 
with  reasonable  accuracy  from  existing  theory.  An  understanding  of  disturbed  orbital 
behavior  permits  possible  control  of  the  elements  and  exploitation  of  the  inevitable 
physical  effects  of  perturbations.  - 

The  magnitude  of  the  effects  of  perturbations  depends  significantly  upon 
the  initial  values  of  certain  orbital  parameters  such  as  semimajor  axis,  eccentricity, 
inclination,  argument  of  perigee,  and  longitude.  Achievement  of  target  parameters 
(passive  control)  can  minimize  deviations  from  the  desired  orbit.  If  needed,  the 
employment  of  onboard  propulsion  resources  (active  control)  may  be  implemented. 
The  consequence  will  be  near  constant  or  beneficially  varying  orbital  parameters — a 


This  thesis  provides  physical  understanding,  theories,  and  mathematical 
formulations  for  different  types  of  frozen  orbits.  Included  in  the  frozen  orbit  concept 
are  minimum  altitude  variation  arcs,  various  categories  of  low-Earth  orbits  such  as 
Sun  synchronous,  and  geosynchronous  orbits.  This  study  serves  as  a  readily 
accessible  and  useful  source  of  information  on  Earth-orbiting  satellites  which  have 
near  constant  or  beneficially  varying  orbital  parameters. 
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INTRODL  CTION 

Since  October  4,  1957,  when  the  Soviet  Union  launched  the  world’s  first 
artificial  satellite,  Sputnik  1,  man-made  objects  have  been  rotating  around  the  Earth  in 
ever-increasing  numbers.  Long  before  this  first  launch,  it  was  known  that  the  actual 
path  that  the  satellite  would  follow  would  be  other  than  pure  Keplerian  motion. 

Theories  endeavoring  to  describe  perturbed  motion  have  been  applied  since 
the  time  of  Newton.  By  his  analysis,  most  of  the  variations  in  the  Moon’s  orbit  were 
explained  in  the  “Principia”  published  in  1687.  Newton  showed  that  the  observed 
inequalities  in  the  motion  of  the  Moon  were  due  to  perturbations  by  the  Sun.  In  1749, 
Clairaut  explained  the  movement  of  the  perigee  of  the  Moon  by  using  second-order 
approximations  [1,9].  A  partial  exposition  of  the  method  of  variations  of  elements  was 
published  by  Euler  four  years  later.  In  1772,  his  perturbation  theory  applied  equations 
of  motion  referenced  to  axes  rotating  with  the  mean  motion  of  the  Moon  [1,6,9]. 

After  publishing  his  technique  of  initial  orbit  determination  in  1780, 
Laplace  turned  his  efforts  towards  the  study  of  perturbed  motion.  His  theory 
transformed  the  equations  of  motion  so  that  the  true  longitude  was  the  independent 
variable.  His  reference  orbit  was  a  Keplerian  ellipse  modified  to  avoid  terms 
proportional  to  the  time. 

In  1783,  Lagrange  determined  that  Newton’s  second-order  system  of 
equations  could  be  written  in  terms  of  an  equivalent  first-order  system.  Generalized 


coordinates  and  Lagrange’s  equations  provide  a  systematic  procedure  for  solving 
systems.  The  resulting  equations  are  equivalent  to  the  equations  of  motion  which 
would  have  been  obtained  by  a  direct  application  of  Newton’s  laws.  However,  the 
Lagrangian  function  contains  only  velocities  and  displacements,  and  no  accelerations 
are  required.  The  equations  produce  a  more  convenient  and  manipulate  form  than  a 
Newtonian  formulation.  When  the  total  perturbing  acceleration  can  be  represented  as 
the  gradient  of  a  scalar  potential  (velocity-independent),  the  resulting  form  of  the 
perturbation  equations  is  known  as  Lagrange’s  planetary  equations  [3, 4, 8, 9],  The 
existence  of  the  perturbing  force  represented  by  the  scalar  potential  causes  the  motion 
to  depart  from  the  unperturbed  motion  and  the  orbital  elements  to  change. 

Gauss,  in  1814,  constructed  an  appropriate  generalization  for  an  arbitrary 
perturbing  force.  It  is  not  always  true  that  the  perturbing  forces  can  be  written  as  the 
gradient  of  a  velocity-independent  scalar.  Atmospheric  drag  is  an  example  of  a 
non-conservative  force.  The  method  resolves  the  perturbing  force  into  three  mutually 
perpendicular  components:  a  component  in  the  direction  of  the  radius  vector,  a 
component  perpendicular  to  the  radius  vector  in  the  orbital  plane,  and  a  component 
perpendicular  to  the  orbital  plane  [1,3,8]. 

After  twenty  years  of  laboring,  Delaunay,  in  1860,  completed  the  most 
exhaustive  application  of  the  canonic  method  in  celestial  mechanics.  The  theory 
involved  the  analytical  removal  of  the  terms  of  the  disturbing  function,  one  by  one. 
Most  working  procedures  based  on  Hamiltonian  dynamics  endeavor  to  eliminate  entire 
classes  of  variables  from  the  original  system  of  differential  equations.While  generating 
suitable  transformations,  the  form  of  Hamilton’s  equations  must  be  preserved  at  all 
times.  This  imposes  constraints  on  the  selection  of  variables.  These  transformations 
are  called  canonical.  The  appearance  of  the  mean  anomaly  in  every  one  of  the  classical 
sets  of  canonic  variables  forces  series  expansions,  not  only  to  accommodate  powers  of 


the  perturbation  parameter,  but  also  to  allow  for  the  transcendental  nature  of  the 
coordinate-time  relation.  The  ability  to  divide  the  Hamiltonian  and  treat  each  term 
separately,  by  a  succession  of  canonic  transformations,  allowed  Delaunay  to  complete 
his  undertaking.  The  solution  extended  to  include  all  terms  of  the  seventh  order  and 
some  of  the  eighth  [2, 4, 5,8]. 

Special  perturbations  refers  to  the  full  numerical  solution  of  the  appropriate 
differential  equations.  One  such  method  was  developed  in  1857  by  Encke.  In  this 
approach,  the  difference  between  the  primary  acceleration  and  all  perturbing 
accelerations  is  integrated.  A  reference  (or  osculating)  orbit  is  established  at  an  epoch. 
The  osculating  orbit  is  the  orbit  that  would  result  if  all  perturbing  forces  were  instantly 
removed.  The  body  at  this  instant  has  the  same  coordinates  and  the  same  velocity 
components  as  in  the  unperturbed  orbit.  At  increments  of  time  later  the  osculating  orbit 
is  compared  to  the  actual  orbit.  When  the  difference  between  the  coordinates  of  the  two 
orbits  becomes  excessive,  a  new  osculating  orbit  is  created  by  a  process  known  as 
rectification.  A  new  epoch  and  starting  point  which  coincides  with  the  true  orbital  path 
is  selected.  A  new  osculating  orbit  is  then  calculated  from  the  true  position  and  velocity 
vector,  and  the  process  of  comparison  is  repeated.  The  advantage  of  Encke’s  method 
is  that  the  difference  between  the  position  vector  for  the  disturbed  motion  and  the 
undisturbed  motion  and  its  derivatives  is  normally  small.  Consequently,  a  large 
integration  interval  is  possible  [1,6,7,10]. 

In  the  early  1900’ s,  P.  H.  Cowell  developed  a  perturbation  method  he 
used  to  determine  the  orbit  of  Jupiter's  eighth  satellite.  The  differential  equations  of 
motion  contain  all  the  forces  acting  on  the  satellite.  The  acceleration  components  are 
integrated  step  by  step  to  produce  velocity  and  position  components.  These  values  are 
then  used  with  the  equations  of  motion  to  calculate  the  acceleration  components.  The 
main  advantage  of  Cowell's  method  is  the  simplicity  of  formulation  and 
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implementation.  Any  number  of  perturbations  can  be  handled  at  the  same  time.  The 
main  disadvantage  of  this  method  is  that  as  the  velocity  increases  near  perigee,  smaller 
integration  steps  must  be  taken.  This  increases  the  number  of  calculations  and 
accumulative  error  due  to  roundoff  [3,6,7]. 

The  mentioned  forerunners  provided  the  necessary  foundation  for  modem 
orbital  mechanics.  By  1957,  the  necessary  theories  were  available  to  amplify  and  to 
study  and  predict  the  effects  of  perturbations. 

The  dominant  force  acting  on  an  Earth-orbiting  satellite  is  central  gravity 
attraction.  Depending  upon  the  altitude,  other  forces  influence  the  motion  in  varying 
degrees.  The  principal  perturbations  result  from  the  Earth’s  oblateness,  the 
longitudinal  variations  in  the  gravity  field,  atmospheric  drag,  lunar-solar  gravitational 
attractions,  and  solar  radiation.  Minor  perturbations  are  also  caused  by  electromagnetic 
interaction  of  a  charged  satellite  with  the  Earth’s  magnetic  field  and  radiation  belts, 
tidal  influences,  re-radiation  of  sunlight  from  the  Earth,  and  relativistic  effects. 
Although  the  perturbations  depend  on  the  direction  and  time  intervals  during  which 
they  act,  the  predominant  disturbing  force  is  due  to  the  Earth’s  oblateness.  As  the 
geocentric  distance  increases,  the  oblateness  effect  diminishes,  the  lunar-solar 
attraction  increases,  the  drag  decreases,  and  the  solar  radiation  pressure  remains 
relatively  constant  [11,12,13]. 

1.1  Related  Literatures 

Applying  the  divergence  theorem  to  Gauss’  Law  and  using  the  fact  that  the 
Newtonian  gravity  field  is  conservative  Poisson’s  equation  is  derived.  At  external 
points  to  the  main  body,  where  the  density  function  is  zero,  Laplace’s  equation  is 
obtained.  The  gravitational  potential  satisfying  Laplace’s  equation  can  be  expressed  in 
terms  of  Legendre's  associated  polynomials.  The  external  gravity  field  can  be  obtained 


from  the  solution  by  differentiation  [4,9,10,12,14]. 

The  gravitational  field  of  the  Earth  can  be  expressed  as  the  sum  of  a  series 
of  spherical  harmonics.  These  are  of  three  types:  (l)the  zonal  harmonics  which 
depend  on  latitude  only,  (2)the  sectorial  harmonics,  which  depend  on  longitude  only, 
and  (3)  the  tesseral  harmonics  which  depend  upon  both  longitude  and  latitude.  Since 
the  Earth  is  rotating  the  tesseral  and  sectorial  harmonics  imply  a  gravitational  field 
which  is  time  dependent  This  is  manifested  as  a  change  in  the  mean  motion. 

The  artificial  satellites  that  followed  Sputnik  1  provided  the  means  to 
investigate  the  Earth’s  gravitational  field.  With  the  data  that  is  obtained  on  the  orbital 
elements  of  satellites,  the  harmonics  are  obtained. 

O’Keefe,  et  al.  used  data  obtained  from  Vanguard  I  (1958)  to  investigate 
the  lower  harmonics.  Making  the  assumption  that  perigee  distance  is  unaffected  by 
drag  and  calculating  the  effects  ol  time-integrals  of  the  drag  corrections  to  the 
semimajor  axis  and  the  eccentricity,  the  effects  of  drag  were  removed  from  the  results. 
It  was  shown  for  this  high  satellite  (perigee  height  =  658  km  and  apogee  height  =  3960 
km)  that  the  orbital  elements  at  37  different  times  over  a  297  day  period  were 
accurately  derived  from  the  zonal  harmonics  of  degrees  0,  2,  3,  and  4  together  with  the 
six  initial  elements  [15]. 

Using  equations  conceived  by  Tisserand  relating  the  ratio  of  semimajor 
axis  and  radius  to  an  integral  of  mean  anomaly,  Kozai  (1957)  derived  the  periodic 
perturbations  of  the  first  order  and  secular  perturbations  up  to  the  second  order  using 
an  averaging  technique  and  mean  elements.  He  assumed  no  air-resistance,  a 
symmetrical  density-distribution  with  respect  to  the  axis  of  rotation,  that  the  coefficient 
of  the  second  harmonic  of  the  potential  is  a  small  quantity  of  the  first  order,  and  those 
of  the  third  and  fourth  harmonics  are  of  the  second  order  [16,17].  Study  of  the  motion 
of  the  Explorer  7  and  third  Vanguard  satellites  produced  coefficients  of  the  second. 


third,  and  fourth  harmonics.  Differences  between  the  observed  and  the  computed  value 
of  the  orbital  elements  implied  the  importance  of  the  higher  harmonics  [18]. 

Care  must  be  taken  in  defining  the  constants  of  integration  when  handling 
the  Lagrange  equations.  During  the  period  1957-1962,  a  large  number  of  theoretical 
papers  were  conceived  to  solve  the  problem  analytically.  The  obvious  choices  are  the 
values  of  the  Keplerian  elements  at  a  particular  instant  of  time.  It  was  convenient  to 
define  the  constants  as  elements  of  a  fictitious  reference  or  intermediate  orbit.  In  some 
theories  this  reference  orbit  is  defined  geometrically,  while  other  theories  defined  the 
constants  dynamically  as  corresponding  to  a  constant  part  of  the  potential  function. 
The  more  familiar  theories  include:  Musen  (1959)  adapted  the  Hansen  lunar  theory  to 
a  form  suitable  for  solution  by  iteration,  Brouwer  (1959)  applied  Von  Zeipel’s  method 
of  canonical  transformation  with  a  purely  Keplerian  intermediary,  Vinti  (1959,1961) 
separated  the  equations  of  motion  by  using  ellipsoidal  coordinates,  King-Hele  (1958) 
employed  a  Keplerian  ellipse  of  fixed  inclination  and  perigee  argument  as  intermediary 
and  solved  in  successive  approximations  according  to  powers  of  the  second  harmonic 
and  eccentricity,  and  Merson  (1961)  developed  a  method  in  a  similar  manner  to 
King-Hele  but  started  from  osculating  elements  when  the  satellite  is  at  the  node  [23]. 
Merson  then  integrated  the  equations  for  the  variation  of  the  osculating  elements  to 
yield  the  complete  perturbations  of  the  first  order  due  to  the  second  harmonic.  Also,  he 
obtained  the  secular  perturbations  of  the  second  order  due  to  the  second  harmonic  and 
of  the  first  order  due  to  the  third  to  sixth  harmonics.  A  set  of  smoothed  elements  was 
derived  in  which  the  perturbations  of  the  even  harmonics  have  no  singularities,  the 
semimajor  axis  and  eccentricity  have  no  variation  due  to  the  second  harmonic,  and  the 
other  elements  have  the  smallest  possible  amplitudes  of  oscillation  [50]. 

Analytical  investigations  of  the  oblateness  effects  showed  the  argument  of 
perigee,  the  ascending  node,  and  the  mean  anomaly  experience  secular  variations  or 


continously  increasing/decreasing  changes  from  the  adopted  epoch  value,  as  well  as 
periodic  variations.  Semimajor  axis,  eccentricity,  and  inclination  experience  only 
periodic  variations.  Periodic  includes  both  short  and  long  periodic  variations  [19,24], 

King-Hele  showed  the  odd  harmonics  determine  the  periodic  changes, 
while  the  even  harmonics  give  rise  to  the  secular  changes.  Data  from  Sputnik  2, 
Vanguard  1,  and  Explorer  7  was  used  to  determine  values  for  the  second,  fourth,  and 
sixth  harmonics  [19].  In  1963,  King-Hele,  et  al.  determined  the  first  seven  even  zonal 
harmonics  in  the  Earth’s  gravitational  potential  using  orbital  information  from  seven 
satellites  covering  a  wide  range  of  latitude  [20].  Transit  2A  and  Samos  2,  two 
high-inclination  satellites,  permitted  a  study  of  inclinations  uniformly  distributed 
between  23°  and  96°.  King-Hele  published  the  updated  values  for  the  even  zonal 
harmonics  in  1964  to  be:  J2  =  (1082.64  ±  .02)  x  10‘6;  J4  =  (-1.52  ±  .03)  x  10  6; 
J6  =  (.57  ±  .07)  x  10  6;  and  Jg  =  (.44  ±  .11)  x  10‘6  [21].  Coefficients  of  the  odd 
zonal  harmonics  were  evaluated  by  analyzing  the  oscillations  in  orbital  eccentricity  of 
fourteen  satellites  covering  the  widest  and  most  uniform  distribution  in  inclination  and 
semimajor  axis  possible  [22]. 

Coefficients  of  some  low  order  tesseral  harmonics  in  the  Earth’s 
gravitational  field  were  successfully  evaluated  in  the  1960s.  The  first  satisfactory 
comprehensive  model  was  the  Smithsonian  Standard  Earth  II  published  in  1970  by  the 
Smithsonian  Astrophysical  Observatory,  Cambridge,  Massachusetts.  This  model 
relied  largely  on  100,000  optical  observations  of  satellites  from  Baker-Nunn  cameras. 
The  expansion  of  the  geopotential  was  truncated  at  degree  and  order  16  which 
produced  about  250  geopotential  coefficients  to  evaluate.  The  geopotential  coordinates 
and  station  coordinates  were  then  adjusted  so  that  the  observations  of  21  satellites  from 
about  30  ground  stations  achieved  the  best  possible  fit  to  these  perturbed  orbits  while 
satisfying  geometrical  constraints  for  simultaneous  observations.  Approximately 
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200,000  equations  were  solved  by  least  squares  [25]. 

In  recent  years  the  Earth’s  gravitational  field  has  been  determined  with 
continually  improving  accuracy,  using  hundreds  of  thousands  of  observations  of  Earth 
satellites,  chiefly  optical,  laser  and  Doppler,  together  with  surface  gravimetry.  More 
recently,  altimeter  measurements  from  satellites  such  as  Geo3  have  been  included  [25]. 

The  launching  of  the  first  artificial  satellite  demonstrated  the  need  to  be  able 
to  estimate  lifetime  of  a  satellite  and  determine  air  density.  Consequently,  the  studies  of 
atmospheric  drag  emerged. 

Equations  expressing  the  variation  of  orbital  period,  perigee  distance,  and 
eccentricity  were  first  developed  by  T.  Nonweiler  (1958)  and  by  King-Hele  (1958). 
Later  King-Hele  extended  his  efforts  to  include  the  effects  of  atmospheric  rotation 
with  oblateness  and  a  realistic  variation  of  air  density  with  height  [29],  The  problem 
was  treated  by  Brouwer  and  Hori  (1961)  on  the  basis  of  a  non-rotating  spherical 
exponential  model  atmosphere.  The  exponential  atmosphere  function  was  expanded 
into  a  series  to  the  fifth  power  which  limited  the  eccentricity  to  a  very  small  quantity. 
For  large  eccentricity,  a  constant  coefficient  for  each  term  of  the  series  was  introduced 
so  that  the  modified  series  fit  the  exponential  function.  The  Brouwer-Hori  theory, 
which  combined  gravitational  and  drag  effects,  was  limited  to  circular  orbits  [26]. 

Instead  of  using  an  exponential  function  for  the  atmospheric  density 
model,  Lane  (1965)  assumed  a  non-rotating  power  function  atmosphere  model. 
Although  Lane  found  it  necessary  to  expand  the  velocity  in  a  series  in  order  to  perform 
the  necessary  integrations,  the  expansion  of  the  density  function  was  not  required. 
Lane  then  used  the  method  of  successive  approximations  to  complete  the  integration. 
His  solution  excluded  orbits  of  small  eccentricity  and  low  inclination  [27].  Later,  Lane 
and  Cranford  (1969)  used  numerical  observations  to  eliminate  the  singularity  of  zero 
eccentricity  and  zero  inclination.  The  latter  three  works  used  the  Von  Zeipel  method  to 


9 

obtain  the  solutions  by  dealing  with  Delaunay  canonical  variables. 

Using  an  asymptotic  method,  a  first-order  approximation  to  the  problem 
was  derived  by  Zee  (1971)  by  dealing  directly  with  the  equations  of  motion  in 
spherical  coordinates  rather  than  with  the  Delaunay  canonical  variables  used  in 
previous  investigations  [30]. 

Two  classifications  of  mathematical  solution  techniques  which  can  be  used 
in  the  study  of  drag  are  numerical  and  analytical.  A  numerical  method  applies 
numerical  integration  to  the  osculating  differential  equations  to  obtain  the  state  at  a  later 
time.  An  analytical  method  generally  uses  analytical  equations  to  transform  from  the 
osculating  to  a  mean  set  of  differential  equations.  The  equations  are  integrated 
analytically  to  predict  a  future  mean  state. 

Recently,  Hoots  (1981)  used  canonical  transformations  and  the  method  of 
averaging  to  obtain  transformations  of  variables.  This  significantly  simplified  the 
transformed  differential  equations.  Analytical  integration  yields  the  six  osculating 
orbital  elements.  The  analytical  solution  is  for  the  motion  of  an  artificial  Earth  satellite 
under  the  influence  of  the  gravitational  zonal  harmonics  through  J4  and  any  dynamic 
atmosphere  [31]. 

The  importance  of  determining  the  lunar  and  solar  perturbations  in  the 
motion  of  an  artificial  satellite  was  demonstrated  by  Kozai’s  discovery  that  certain 
long-period  terms  in  the  development  of  the  disturbing  function  cause  large 
perturbations  in  the  orbital  elements.  These  perturbations  may  strongly  affect  the 
perigee  height  and  the  lifetime  of  a  satellite.  Modifying  A.  Cayley’s  development  of  the 
solar  perturbative  function,  Musen  et  al.  generated  relations  between  perigee  height 
variations  and  launch  conditions.  Using  data  from  Vanguard  I  and  Explorer  VI  values 
of  perturbations  in  the  perigee  height  were  computed  using  a  trigonometric  expansion 
of  the  disturbing  function  with  the  angular  equatorial  elements  of  the  Sun  and  Moon  as 


arguments.  The  disturbing  function  is  then  integrated  by  approximating  the  equatorial 
elements  as  linear  functions  of  time.  Kaula  adapted  the  mathematical  development  of 
Musen  to  a  form  convenient  for  use  in  connection  with  the  analysis  of  satellite  orbits 
affected  by  a  terrestrial  gravitational  field.  Kaula  further  expressed  the  disturbing 
function  in  terms  of  osculating  Keplerian  elements.  The  formulae  use  equatorial 
elements  for  the  Moon.  The  disturbing  function  yields  all  significant  first-order 
lunar-solar  effects  when  used  with  the  equations  of  motion  [12,33].  Neither 
inclination  nor  node  of  the  Moon’s  orbit  with  respect  to  the  equator  of  the  Earth  are 
simple  functions  of  time.  The  same  elements  with  respect  to  the  ecliptic  are  nicely 
approximated  by  a  constant  and  a  linear  function  of  time  respectively.  Giacaglia  ( 1 973) 
presented  a  paper  in  which  he  obtained  the  disturbing  function  for  the  lunar 
perturbations  using  ecliptic  elements  for  the  Moon  and  equatorial  elements  for  the 
satellite.  The  secular,  long-period,  and  short-period  perturbations  are  then  computed 
in  closed  form  in  both  inclination  and  eccentricity  of  the  satellite  [34], 

The  motion  of  a  satellite  under  the  influence  of  the  longitude  dependent 
terms  of  the  geopotential  and  in-plane  station-keeping  requirements  of 
geosynchronous  satellites  with  nearly  zero  eccentricity  and  inclination  were 
investigated  by  Blitzer  et  al.  (1962,196-.  Izsak  (1961),  Frick  (1962),  Allan  (1963), 
Wagner  (1964,1965),  Blitzer  (1965),  and  others  [39-46].  Their  findings  regarding 
equilibrium  positions  and  librationa!  periods  are  in  harmony. 

The  study  of  a  satellite  with  a  mean  motion  commensurable  with  the 
Earth’s  rotation  was  generalized  to  include  the  effects  of  tesseral  harmonics  on 
near-circular  orbits  of  arbitrary  inclinations  in  the  works  of  Blitzer  (1966),  Cook 
(1966),  and  Wagner  (1966)  [47-49]. 

Ordinarily,  the  solar  radiation  pressure  has  a  very  small  effect  on  the  orbit. 
This  perturbation  becomes  more  important  for  geosynchronous  satellites  with  large 


solar  arrays.  The  major  effects  are  the  oscillation  of  the  eccentricity  vector  and  the 
variation  in  perigee  height. 

Kozai  (1961)  developed  closed  form  analytical  expressions  to  include 
short-term  periodics  for  first  order  perturbations.  The  three  components  of  the 
disturbing  equations  were  based  on  three  assumptions:  (1)  the  parallax  of  the  Sun  is 
negligible,  (2)  the  solar  flux  is  constant  along  the  orbit  of  the  satellite  if  there  is  no 
shadow,  and  (3)  there  is  no  re-radiation  from  the  surface  of  the  Earth.  The  resulting 
expressions  ire  then  evaluated  over  two  limits  of  the  independent  variable,  eccentric 
anomaly,  which  must  be  derived  by  numerical  methods  [35], 

The  short-term  secular  variations  in  period  resulting  from  solar  radiation 
pressure  were  studied  by  Wyatt  (1961).  It  was  shown  that  the  effect  of  solar  radiation 
is  negligible  when  the  satellite  is  in  continuous  sunshine.  During  the  time  the  satellite  is 
in  the  Earth’s  shadow,  the  secular  acceleration  may  attain  substantial  values,  positive 
or  negative  depending  on  the  orientation  of  the  orbit  relative  to  the  Sun.  Wyatt 
presented  a  general  formula  for  computing  secular  accelerations  as  far  as  terms  in  the 
square  of  the  eccentricity  [36J. 

Based  on  the  assumption  that  the  disturbing  acceleration  is  constant  while 
the  satellite  is  in  the  sunlight  and  is  zero  in  the  Earth’s  shadow.  Cook  (1962)  generated 
equations  for  the  variations  in  orbital  elements.  The  perturbation  expressions  are 
formulated  in  terms  of  the  radial,  transverse,  and  normal  components  of  the  disturbing 
force  evaluated  at  perigee  [38]. 

An  expression  which  favorably  compares  to  observations  was  presented 
by  Soop  for  the  mean  drift  rate  of  the  eccentricity  vector  as  a  function  of  radiation 
pressure,  velocity,  mass,  cross-sectional  area,  and  sidereal  angle  of  the  Sun  [37], 


Over  the  last  20-25  years,  theories  to  treat  the  major  perturbations  have 
been  developed  and  used  with  much  success.  This  study  has  endeavored  to  generalize 
the  frozen  orbit  concept  to  include  types  of  orbits  that  have  near  constant  or 
beneficially  varying  orbital  parameters.  To  achieve  these  frozen  orbits  the  effects  of 
perturbations  are  exploited  to  maintain  the  type  of  orbit  desired. 

Chapter  II  presents  the  equations  of  motion  for  a  perturbed  orbit,  the 
Lagrange  Planetary  Perturbation  Equations  explicitly  expressed  in  terms  of  the 
disturbing  forces,  and  analytical  expressions  for  the  perturbations  due  to  an  axially 
symmetric  gravitational  field.  Secular,  long-term  periodic,  and  short-term  periodic 
variations  are  distinguished.  Secular  variations  of  the  second  order  are  expanded  and  a 
root  finder  is  used  on  the  resulting  high  order  polynomials  to  find  the  eccentricity  and 
inclination  that  freezes  the  argument  of  perigee  and  the  ascending  node  separately.  The 
orbital  elements  which  freeze  the  perigee  along  the  meridian  are  studied.  Using  Bessel 
expansions  the  change  in  radius  is  investigated  at  any  point  on  the  orbit  including  the 
perigee.  The  chapter  concludes  with  a  discussion  of  a  frozen  orbit  where  the  change  in 
eccentricity  due  to  the  third  harmonic  is  zero  as  is  the  change  in  argument  of  perigee 
due  to  the  second  and  third  harmonics. 

Atmospheric  drag  effects  are  investigated  in  Chapter  III.  The  force  which 
drains  energy  from  the  satellite  and  causes  the  orbit  to  shrink  is  analyzed  for  orbits  of 
normal  eccentricity  (0  <  e  <  0.2).  Changes  in  orbital  elements  are  handled  both 
analytically  and  numerically  usi.  _  a  fourth-order  Runge-Kutta  scheme.  Scale  height  is 
determined  using  density  values  from  U.S.  Standard  Atmosphere  tables.  Between 
cardinal  altitudes  an  exponential  atmosphere  is  assumed. 

Minimum  altitude  variation  arcs  are  dealt  with  in  Chapter  IV.  A  least 
squares  method  is  derived  for  determining  orbital  elements  for  an  orbit  that  closely 


follows  the  contour  of  the  oblate  Earth  over  fixed  latitude  ranges.  Results  are  included 
for  both  drag  and  no-drag  effects.  Impulse  requirements  to  maintain  the  orbit  are 
treated. 


Sun  synchronism  is  essential  to  maintain  the  same  lighting  conditions  over 
a  particular  area  of  the  Earth.  Orbital  parameters  that  are  necessary  to  achieve  this  type 
of  frozen  orbit  are  specified  in  Chapter  V.  The  chapter  considers  the  ground  trace  of  an 
artificial  satellite.  An  analytical  expression  which  includes  oblateness  effects  is 
developed  determining  the  semimajor  axis  that  is  needed  to  give  a  repeated  ground 
trace  in  the  desired  number  of  days. 

The  last  type  of  frozen  orbit  analyzed  is  the  geosynchronous  orbit.  Chapter 
VI  includes  the  main  perturbing  effects  on  this  type  of  orbit  and  longitude  station 
keeping  requirements.  The  chapter  concludes  with  a  brief  discussion  of  the  effects  of 
solar  radiation. 

1.3  Contribution 

This  thesis  provides  physical  understanding,  theories  and  mathematical 
formulations  for  different  types  of  frozen  orbits.  It  serves  as  a  convenient  and  useful 
source  of  information  on  Earth  orbiting  satellites  which  have  near  constant  or 
beneficially  varying  orbital  parameters. 


CHAPTER  n 


GRAVITATIONAL  PERTURBATIONS  ON  ORBITAL  PARAMETERS 


2.1  Equations  of  Motion 

The  equations  of  motion  for  the  classical  two-body  problem  are  governed 
by  the  inverse  square  law.  The  orbit  may  be  determined  from  the  equations  of  motion 
and  initial  values  of  position  and  velocity.  In  the  two-body  problem,  all  forces  are 
neglected  except  the  mutual  gravitational  attraction  of  two  spherically  symmetrical 
bodies.  The  resulting  basic  differential  equation  is 


(2.1) 


This  is  equivalent  to  three  second-order  differential  equations  requiring  six 
constants  of  integration  for  the  complete  solution.  The  solution  may  be  obtained  by  a 
variety  of  methods.  The  result  is  the  general  polar  equation  for  a  conic  section  with  the 
origin  at  a  focus  [1,51]: 


1  +  e  cos  v 


(2.2) 


Equation  (2.2)  when  combined  with  the  constant  angular  momentum  can  provide  the 
last  integration,  the  time  relation. 
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In  actuality,  forces,  in  addition  to  the  central  gravitational  force,  act  on  the 
system.  Other  forces  include  Earth  oblateness,  atmospheric  resistance,  lunar-solar 


attractions,  solar  radiation,  etc.  The  equation  of  motion  becomes 


+  f 


(2.3) 


In  some  cases,  it  will  be  convenient  to  work  with  the  total  perturbing 
force.  However,  if  the  force  is  conservative,  it  can  be  represented  as  the  gradient  of  a 
scalar  function.  Equation  (2.3)  becomes 

r  +  JLL  =VR  (2-4) 

r3 

where  R  is  the  disturbing  function.  The  form  of  R  will  depend  on  the 
particular  type  of  perturbing  source.  For  multiple  perturbing  sources,  the  respective 
disturbing  functions  can  be  added  linearly  to  yield  the  total  R  [2,4,8,10). 


2.2  Lagrange’s  Planetary  Equations 

If  a  set  of  orbital  elements  (a,  e,  i,  Q,  co,  M)  is  obtained  at  a  fixed  instant, 
an  osculating  orbit  is  defined.  The  osculating  ellipse  is  the  path  that  would  be  followed 
by  the  satellite  if  the  perturbing  forces  were  all  suddenly  removed.  The  satellite  at  this 
time  has  the  same  coordinates  and  the  same  velocity  components  in  the  perturbed  as  in 
the  unperturbed  orbit.  Stated  differently,  the  satellite  has  the  same  position  and  velocity 
as  it  would  in  purely  two-body  motion.  If  an  osculating  ellipse  corresponding  to  a 
later  time  is  considered,  the  differences  in  the  orbital  elements  from  those  of  a  previous 
time  are  caused  by  the  perturbations.  One  advantage  in  expressing  actual  orbital  motion 


in  terms  of  osculating  elements  is  that  although  the  position  and  velocity  vectors 
change  rapidly  with  respect  to  time,  the  Keplerian  elements  vary  slowly  even  in 
perturbed  motion. 

The  Lagrange  planetary  equations  provide  the  means  for  finding  the  rates 
of  change  of  the  osculating  elements.  Since  any  six  linearly  independent  combination 
of  orbital  elements  provides  a  valid  set,  the  equations  may  be  written  in  a  variety  of 
forms.  For  this  undertaking  the  perturbing  force  f  will  be  expressed  in  terms  of  the 
following  three  components:  (1)  fj  is  along  r,  (2)  f2  is  perpendicular  to  r  in  the 
osculating  plane  and  in  the  direction  of  increasing  anomaly,  and  (3)  f3  is  perpendicular 
to  the  osculating  plane  in  the  direction  of  the  angular  momentum  vector.  Substituting 
angular  momentum  in  the  equivalent  terms,  a  modified  form  of  Lagrange's  planetary 
equations  is  [8,50]: 


dT  =2(ir)  Vsmv+fji) 


df  =  (jrHfi  sinv  +  f2cosv  (1  +-£-)  +  f2^-l 


r  f3  cos  u 
h 


dfl  _  £ h  sin  u 
dt  -  (h  sin  i) 


df =  (pe  ^'fi  cosv  +  f2sin  v  (l+_p)-f3e  sin  u  cot  i  (-£-)] 


df  =  W  f  l (cos  v  -  ^§-r  )  -f2  sin  v  (i+X)] 


(2.10) 


'/£y«y*\sV*V  ,,w’* 


LIST  OF  SYMBOLS 


a  =  semimajor  axis 

e  =  eccentricity 

E  =  eccentric  anomaly 

f  =  acceleration  vector  of 
perturbing  force 

fb  f^  f3=  components  of  f 

F  =  perturbing  potential 

h  =  (jip)1'/2=  angular  momentum 
for  an  ellipse 

i  =  inclination  angle 

J  n  =  zonal  harmonics 

k  =  sin  2  i 

M  =  mean  anomaly  =  nt  +  cr 
n  =  mean  motion  =  (^/a3 
p  =  latus  rectum  =  a  (1  -  e2 ) 


Pn  =  Legendre  polynomial  of  order  n 

r  =  radius  vector 

r  =  magnitude  of  r 

R  =  mean  equatorial  radius  of 
the  Earth 

U  =  gravitational  potential 
u  =  v  +  0) 


p  =  gravitational  constant  = 


Q  =  right  ascension  of  the  node 


v  =  true  anomaly 
0  =  colatitude 


co  =  argument  of  perigee 
a  =  mean  anomaly  at  the  node 
a  =  modified  mean  anomaly  = 

M - f  1  ndt 
J  o 
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Figure  2.1  ORBITAL  PARAMETERS 
p  =  perigee,  s  =  satellite  position 

2.3  ^.Gravitational  Field 

The  force  exerted  by  a  conservative  force  field  acting  on  a  particle  (satellite) 
can  be  written  as  the  gradient  of  gravitational  potential  U.  In  terms  of  equation  (2.4), 
VU  =  VR. 

If  variations  with  longitude  are  ignored  and  the  Earth  is  assumed  to  be 
axially  symmetric,  the  gravitational  potential  U  at  an  exterior  point  distance  r  from  the 
center  of  the  Earth  and  with  a  colatitude  0  may  be  written  as  a  series  of  spherical 
harmonics  in  the  form  [4,  9,  1 1,  12,  19] 

U  =  (JtL)  [  1- 1  Jn  (4)n  Pn(cos0)] 

r  n=2  '  n 


(2.11) 
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The  principal  perturbation  on  a  near-Earth  satellite  is  due  to  the 
axially-symmetric  (zonal)  harmonics  of  the  Earth’s  gravity  field.  The  perturbing 
potential  F  for  this  perturbation  becomes: 

F  =  U-T  =  -  Or)  nVn(f)npn(cos0>  (2-12) 

The  component  of  the  perturbing  acceleration  f  in  any  direction  is  the  rate 
of  change  of  F  in  that  direction.  Using  spherical  trigonometry  based  on  figure  2.1  it 
follows  that  fj  =  0F/gr 

f2  ,  .(COS  u  Sin  i  /rsjn  e)3F/ae  <2.13) 

fJ  =  -(COSi/rsin8)8F'ae 

Substituting  the  explicit  formulae  for  Pn  up  to  n=6,  performing  the 
differentiations,  then  setting  cos  9  =  sin  i  sin  u,  and  using  the  abbreviations 
S  =  sin  u  and  k  =  sin2i  ,  the  result  is: 


fj  =  ((i/r2)  [3/2  J2  (R/r)2  (3kS2  -  l)+2  J3  (R/r)3sin  i  (5kS3  -  3S) 

+  5/g  J4  (R/r)4  (35k2S4  -  30kS2  +  3) 

+  (3/4)  J5  (R/r)5  sin  i  (63k2S5  -  70kS3+  15S) 

+  ?/i6J6  (R/r)6(231k3S6-315k2S4+  105kS2  -  5)]  (2.14) 

f2  =  -{\\Jr-)  sin  i  cos  u  [3  J2  (R/r)2  sin  i  sin  u  +  (3/2)  J3  (R/r)3  (5kS2  -1)+ 
5/2  J4  (R/r)4  sin  i  (7kS3  -  3S)  + 

'%  J5  (R/r)5  (2lk2S4  -  14kS2  +  1)  + 

21/g  J6  (R/r)6sin  i  (33k2S5  -  30kS3  +  5S)]  (2.15) 


f 3  =  -  (fi/r2)  cos  i  x  [  same  as  f2  ] 


(2.16) 
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To  determine  the  secular  and  long-period  terms  equations  (2.14),  (2.15), 
and  (2.16)  are  substituted  in  equations  (2.5)-(2.10).  The  resulting  equations  are  then 
multiplied  by  an  appropriate  equality  with  c^t/£ju  to  obtain  J2  to  J6  and  J22  terms.  This 
allows  the  argument  of  latitude  to  be  used  as  the  independent  variable.  Expressing  the 
derivative  of  the  orbital  elements  with  respect  to  the  argument  of  latitude  as  the  sum  of 
a  number  of  terms  of  a  general  function,  utilizing  appropriate  identities,  and  integrating 
from  one  ascending  node  to  the  next  (u=0  to  u=27t),  the  desired  equations  were 
obtained  (for  details  see  ref.  50,  p  23-31). 

These  analytical  expressions  reveal  that  the  elements  to,  ft,  and  M 
experience  secular  variations  from  the  adopted  epoch  values  as  well  as  periodic 
variations.  The  elements,  a,  i,  and  e  possess  only  periodic  variations  and  oscillate 
about  their  mean  values.  The  long-period  variations  are  caused  by  the  continuous 
variance  of  co.  The  strictly  secular  terms  contain  only  even  Jn  terms  while  the  odd 
zonal  harmonics  contribute  to  periodic  terms.  Since  the  changes  in  orbital  elements 
over  one  period  are  to  be  studied  here,  the  secular  variations  are  of  primary  concern. 
Consequently,  the  following  equations  are  absent  of  any  odd  Jn  terms  since  they  all 
have  products  of  trigonometric  functions  of  co  or  multiples  of  co. 

Aa  =  Ae  =  Ai  =0 

see  sec  sec 


Aftsec  =27t[lJn(R/p)nftn  +  J2(R/p)4ft22]  (2.17) 

ft2  =  -3/2  cos  i 

ft4  =  15/4  cos  i  (  (1  -  7/4  k)  (1  +  3/2e2)  ] 

ft6  =  -105/,6  cos  i  [  (1  -  9/2  k  +  33/g  k2)  (1+  5e2  +  15/g  e4)  ] 

ft22  =  3/2  cos  i  |  3/4  -  5k  -  (*/4  +  5/16  k)  e2] 
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Acosec  =  2 n  [XJ  (R/p)n  &)n  + J.2(R/p)4co„] 


(2.18) 


co2  =  3  (1  -  5/4  k) 

<o4  =  - 15/32  [  (16  -  62  k  +  49  k2)  +  (18  -  63  k  +  189/4  k2)  e2  ] 
w6  =  525/64  [  8/5  (1  -  8  k  +  12%  k2  -  297/32  k3) 

+  6(1-  43/6  k  +  109/g  k2  -  121/g  k3)  e2 
+  (2  -  27/2  k  +  "/4  k2  -  429/32  k3)  e4  ] 
co22  =  (9/4)  [  (95/12  k  -  445/48  k2)  +  (7/12  -  3/gk  -  15/32  k2)  e2  ] 


A  =  n  +  2  n [XJ n  (R/p)n  a„  +  J 2(R7p)4 o*22l 
a* 2  =  3/2(l  -e2)1/2(l  -3/2  k) 


(2.19) 


a  *4  =  45/16  (1  -  e2)1/2  [  (-1  +  5k  -  35/g  k2)  e2  ] 
a  *6  =  -35/i6  (1  -  e2)1/2  [  (1  -  21/2  k  +  189/g  k2  -  231/16  k3) 
x  (1  -  5/2  e2  -  15/g  e4) ) 
a*22  =  9/2(l-e2)-1/2l(25/12k-131/48  k2) 

+  (S/12  ‘49/12k  +  67/48  k2>  e2l 


A  current  procedure  that  is  used  to  simplify  a  system  of  nonlinear 
differential  equations  is  the  method  of  averaging.  By  proper  selection  of  new  variables 
the  original  equations  are  transformed  into  a  simpler  form,  free  of  certain  variables. 
The  normal  form  for  averaging  includes  a  set  of  slow  variables  with  time  variations 
proportional  to  a  small  parameter.  The  fast  variables  have  dominant  parts  with  time 
variations  proportional  to  time.  The  transformation  of  variables  removes  the 
dependence  on  mean  anomaly,  the  fast  variable.  Since  the  method  of  averaging 
requires  perturbations  to  be  continuous  and  periodic  in  the  fast  variable,  it  is  an  ideal 


technique  to  use  to  analyze  orbital  motion.  The  solution  for  these  transformed 
equations  permits  the  calculation  of  the  six  osculating  orbital  elements. 


Expressing  the  disturbing  function  explicitly  in  terms  of  the  orbital 
elements,  the  slow  variable  perturbation  equations  have  the  form: 
n  =  egjCn,  e,  i,  to,  ft,  M) 
e  =  eg2(n,  e,  i,  co,  ft,  M) 

di/dt  =  e>  h  M) 

co  =  eg4(n,  e,  i,  co,  ft,  M) 
ft  =  eg5(n,  e,  i,  to,  ft,  M) 
e  is  a  small  parameter  of  order  J2. 

The  fast  variable  equation  is  of  the  form: 

M  =  n  +  euj(n,  e,  i,  to,  ft,  M) 

gj  thru  g5  and  ut  are  continuous  functions  of  the  orbital  elements  as  well 
as  periodic  functions  of  mean  anomaly  with  a  period  of  2k.  The  equations  are  in  the 
normal  form  for  the  method  of  averaging  and  the  procedure  outlined  in  references  52 
and  54  can  be  employed. 

Using 

F  =  •—  £  K  (~)  Pn(cos9) 

i  n=2  1 

and  adhering  to  the  prescribed  approach,  the  periodics  can  be  eliminated  and  the 
following  secular  variations  of  the  first-order  harmonic  J2  and  the  second-order 


harmonics  J3  and  J4  result: 


(2.20) 


n  =  0  (Since  n  =  -  3n/2a  a,  a  =  0.) 
e  =  0 
di/  -  0 

'dt  ~  u 

0)  =  3/4  n  J2  (R/p)2  (4  -  5k) 

- 15/32  n  J4  (R/p)4  [  16  -  62k  +  49  k2  +  3/4  (24  -  84  k  +  63  k2)  e2  ] 

+  (3/16)n  J22(R/p)4  [72  -169k  +  395/4  k2  +  (-5  +  57/2k  -  225 /g  k2)e2] 

=  -3/2  n  J2  (R/p)2  cos  i  +  15/16  n  J4  (R/p)4  cos  i  [(4  -  7k)  (1  +  3/2  e2)] 

-  3/2  n  J22  (R/p)4  cos  i  [%  +  3/2  (l-e2)1/2  -  k  (5/2  +  %  (1  -  e2)1/2) 

+  e2/4  (1  +  5/4k)  ] 

M  =  n  [  1  +  3/2  J2  (R/p)2  (1  -  3/2  k  (l-e2)1/2)] 

-  45/i28 n  J4  (^P)4  (8  -  40k  +  35k2)  e2  (l-e2)1/2 

+  %  n  J22  (R/p)4  (l-e2)1/2  [  3/2  -  47/12  k  +  139/48  k2 
+  (-V24+13/,2k  +10/96k2)  e2] 

Although  the  J22  term  is  slightly  different,  these  results  corroborate  those 
formulated  by  Merson. 

2.5  The  Argument  of  Perigee 

One  of  the  most  important  perturbations  is  the  rotation  of  the  major  axis  of 
an  orbit  in  its  own  plane.  From  equation  (2.18)  it  is  seen  that  on  the  0(J2)  the  change 
in  argument  of  perigee  over  one  period  is: 

Aco  =  (3rt/2)  J2  R2  (4  -  5  sin2i)  /  (a2(l-e2)2)  (2.21) 

The  angle  at  which  the  argument  of  perigee  freezes  is  normally  referred  to 
as  the  critical  inclination.  For  Aco  =  0,  4  -  5  sin2  i  =  0.  This  occurs  at  an  inclination 
angle  of  63.4'  for  prograde  orbits  and  1 16.6'  for  retrograde  orbits.  For  i  <  63.4°  the 
perigee  moves  in  the  same  direction  as  the  satellite’s  motion;  for  63.4°  <  i  <  1 16.6° 
the  perigee  moves  in  the  opposite  direction;  and  for  i  >  1 16.6°  the  perigee  moves  in 


the  same  direction. 

If  gravitational  perturbations  due  to  the  geopotential  alone  are  considered 
and  if  tesseral  and  sectorial  harmonic  resonances  are  avoided,  only  the  zonal 
harmonics  cause  secular  and  long-period  perigee  fluctuations.  Studies  where  the 
averaged  Hamiltonian  was  expanded  about  critical  inclinations  revealed  very  slow  (on 
the  order  of  a  century)  simple  pendulum-type  oscillations  in  argument  of  perigee  [53]. 

To  improve  the  accuracy  to  0(J22)  the  entire  equation  (2.18)  can  be  set 
equal  to  zero.  Since  there  is  no  secular  change  in  a,  e,  and  i  caused  by  the  Earth’s 
geopotential  these  values  can  be  assumed  constant.  Using  a  binomial  expansion  the 
following  sixth-order  polynomial  equation  in  sin  i  results: 

o  =  [3  J2  a4  (e8  -  4e6  +  6e4  -  4e2  +  1)  (2.22) 

-  ( 1 5  J4  R2  a2  /  32)  ( 1 8  e6  -  20  e4  -  14e2+  16) 

+  524/m  J6R4  (32/20+6e2+2e4)+9/384  J22  a2R2  (56e6-l  12e4+56e2)  ] 

+  [3J2a4/4)(-5e8+20e6-  30e4+  20e2-  5)-,5/32  J4R2a2(-63e6+  64e4+  61e2-  62) 

-  52J'„  R4  (-45  V4  • •  2%> 

+  %84  J22  a2  r2  (*36e6  +  832e4  -  l,556e2  +  760)  ]  sin2i 

+  [  (-15  J4  R2  a2  /  32)  (189/4e6  -  91/2  e4  -  203/4e2  +  49) 

^  524/64  J6r4  (792/32  e4  +  327/4  e2  +  5%) 

+  %84  J22  a2  r2  (‘45e6  -  800e4  +  U35e2  -  890)  ]  sin4i 

-  524/64  J6  R4  (429/32e4  +  363/4e2  +  297/20)  sin6> 

Table  2-1  shows  the  values  of  inclination  that  will  freeze  the  argument  of 
perigee  given  an  eccentricity  and  semimajor  axis.  The  data  reveals  that  as  semimajor 
axis  increases  the  inclination  angle  required  decreases.  The  decrease  is  small  when  the 
eccentricity  is  small.  As  the  eccentricity  increases  the  decrease  in  inclination  is  of 
greater  magnitude  over  the  same  range  of  semimajor  increase.  The  other  trend  that  is 
obvious  is  that  for  a  given  semimajor  axis,  as  the  eccentricity  increases  the  inclination 
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required  decreases. 


A  comparison  of  computed  change  in  argument  of  perigee  with  actual  1984 
data  from  NORAD  confirms  the  contribution  of  zonal  harmonics  of  0(J2)  to  the 
secular  variation  (see  table  2-2).  In  all  cases  the  secular  equation  produces  reasonable 
accuracy.  Deviation  that  occurs  is  possibly  due  to  periodic  variations  caused  by  the 
gravitational  field  as  well  as  changes  due  to  other  major  perturbing  forces.  The  most 
error  per  time  lapse  occurs  with  the  satellite  experiencing  the  most  drag  effect.  Another 
source  of  error  is  the  fact  that  the  change  is  averaged  over  the  entire  time  as  opposed  to 
adding  the  change  that  occurs  each  orbit  over  the  time  lapse. 

Table  2-3  shows  the  first  and  second-order  changes  in  the  argument  of 
perigee  given  a  semimajor  axis,  an  eccentricity,  and  an  inclination  angle.  The  data 
demonstrates  that  for  a  given  eccentricity,  as  the  semimajor  axis  increases  the  rate  of 
change  in  argument  of  perigee  decreases.  Also,  at  a  constant  semimajor  axis,  as  the 
eccentricity  increases  the  rate  of  change  in  perigee  is  insignificantly  altered. 


Table  2-2.  ACTUAL  VS.  COMPUTED  CHANGE  IN  PERIGEE 


Satellite  # 

a 

e 

i 

At 

^actual 

Aco 

compi 

(km) 

(deg) 

(days) 

(deg) 

(deg) 

62  B-P 

6864.903 

.00105 

90.613 

1.90128 

-5.3478 

-7.3246 

62  B-M 

7506.838 

.00702 

50.150 

4.79179 

14.2662 

14.1413 

77-79  F 

7844.892 

.0005 

74.026 

30.42327 

-39.8920 

-45.5713 

19-41  A 

6934.313 

.0024 

55.016 

46.06976 

120.2862 

110.2657 

1 1416  U 

7187.775 

.0012 

98.570 

10.11341 

-28.8684 

-29.4842 

12553 

7227.384 

.00133 

99.026 

4.39037 

-12.9535 

-12.3858 

Table  2-3  CHANGE  IN  ARGUMENT  OF  PERIGEE 


e  =  .02  a  =  7,000  km 


i  (deg) 

0(J2) 

(deg/day) 

0(J22) 

(deg/day) 

diff 

(deg/day) 

30 

9.9013 

9.9318 

-.0304 

35 

8.4793 

8.4939 

-.0146 

40 

6.9637 

6.9669 

-.0031 

45 

5.4007 

5.3969 

.0038 

50 

3.8736 

3.8307 

.0069 

55 

2.3221 

2.3149 

.0071 

60 

.9001 

.8942 

.0058 

65 

-.3851 

-.3891 

.0039 

70 

-1.4946 

-1.4969 

.0023 

75 

-2.3945 

-2.3957 

.0012 

80 

-3.0576 

-3.0583 

.0006 

85 

-3.4637 

-3.4642 

.0004 

90 

-3.5983 

-3.5984 

.0001 

e  =  .02  a  =  7,500  km 


7.777 

6.660 

5.469 

4.242 


7.8088 

6.6734 

5.4688 

4.2317 

2.9993 

1.8084 

.6940 

-.3111 

-1.1774 

-1.8793 

-2.3962 

-2.7126 

-2.8148 


.08  a  =  7,500  km 


Wi 

! 

Bp 

•  Hij|| 

1 

1 

i 

«  PBI 

*  ^Bllt 

1 

Another  strong  secular  perturbation  of  major  significance  is  the  steady 
rotation  of  the  orbital  plane  about  the  Earth’s  axis  commonly  referred  to  as  regression 
of  the  nodes.  Equation  (2.17)  shows  that  all  secular  terms  contain  a  product  of  cos  i. 
Because  of  the  major  influence  of  the  J2  term,  for  i  <  90°,  the  orbital  plane  rotates 
about  the  Earth’s  axis  in  the  direction  opposite  to  the  satellite  motion.  A  90°  inclination 
(polar)  orbit  would  freeze  the  ascending  node  and  there  would  be  no  secular  motion. 

Table  2-4  shows  the  rate  of  change  of  the  ascending  node  as  a  function  of 
eccentricity,  semimajor  axis,  and  inclination  angle.  As  the  eccentricity  increases  for  a 
given  semimajor  axis  and  inclination,  the  rate  of  node  change  increases.  For  a  given 
eccentricity  and  inclination,  an  increase  in  semimajor  axis  causes  a  decrease  in  the  rate 
of  node  change.  Finally,  the  difference  between  first  and  second-order  solutions 
becomes  smaller  when  any  of  the  following  occurs:  eccentricity  decreases,  semimajor 
axis  increases,  or  inclination  increases. 

A  comparison  is  made  between  actual  and  computed  change  in  the 
ascending  node  in  Table  2-5.  The  1984  NORAD  data  is  on  the  same  satellites  listed  in 
Table  2-2. 
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An  interesting  combination  of  the  change  in  argument  of  perigee  and 


ascending  node  is 

£2  +  co  cos  i  =  0  (2.23) 

When  this  relationship  is  satisfied,  the  perigee  is  frozen  along  the  meridian. 
For  all  possible  eccentricities  and  all  semimajor  axes  as  large  as  15,000  km  the 
inclination  (other  than  90°)  for  which  equation  (2.23)  is  true  varies  between  38°  and 
39.2°.  Figure  2.2  shows  the  change  in  ascending  node,  change  in  perigee,  and 
Q  +  co  cos  i  graphically. 


Table  2-4  CHANGE  IN  THE  ASCENDING  NODE 


e  =  .02 

a  =  7,000  km 

i (deg) 

0(J2) 

(deg/day) 

0(J22) 

(deg/day) 

diff 

(deg/day) 

30 

-6.2362 

-6.2515 

.0153 

40 

-5.5162 

-5.5258 

.0095 

50 

-4.6287 

-4.6354 

.0066 

60 

-3.6004 

-3.6063 

.0058 

70 

-2.4628 

-2.4680 

.0051 

80 

-1.2504 

-1.2536 

.0031 

90 

.0000 

.0000 

.0000 

e  =  .05 

a  =  7,000  km 

30 

-6.2625 

-6.2780 

.0155 

40 

-5.5395 

-5.5491 

.0096 

50 

-4.6482 

-4.6549 

.0067 

60 

-3.6156 

-3.6215 

.0058 

70 

-2.4732 

-2.4784 

.0052 

80 

-1.2557 

-1.2589 

.0032 

90 

.0000 

.0000 

.0000 

i (deg) 


a  =  7,500  km 


0(J2) 

(deg/day) 


-4.8983 

-4.3328 

-3.6357 

-2.8280 

-1.9345 

-.9821 

.0000 


-4.9086 

-4.3395 

-3.6404 

-2.8321 

-1.9379 

-.9842 

.0000 


a  =  8,000  km 


(deg/day) 


30 

-4.1407 

-4.1491 

40 

-3.6626 

-3.6681 

50 

-3.0733 

-3.0771 

60 

-2.3906 

-2.3937 

70 

-1.6352 

-1.6379 

e  =  .06 

a  =  12,000  krr 

i 

i (deg) 

0(J2) 

(deg/day) 

0(J22) 

(deg/day) 

diff 

(deg/day) 

20 

-1.5111 

-1.5140 

.0028 

40 

-1.2318 

-1.2332 

.0013 

60 

-.8040 

-.8046 

.0006 

80 

-.2792 

-.2795 

.0002 

e  =  .06  a  =  15,000  km 


20 

-.4728 

-.4730 

.0002 

40 

-.3854 

-.3856 

.0001 

60 

-.2515 

-.2516 

.00009 

e  =  .54 

a  =  15,000  km 

20 

-.9354 

!  -.9371 

.0017 

40 

-.7625 

!  -.7633 

.0007 

60 

-.4977 

|  -.4980 

.0003 

Table  2-5 

ACTUAL  VS.  COMPUTED  CHANGE  IN  ASCENDING  NODE 


Satellite  # 

At 

days 

AQ 

(deg)  actual 

AQ 

(deg)  computed 

62  B-P 

1.40128 

.1593 

.1639 

62  B-M 

4.79178 

-17.3097 

-17.3221 

77-79F 

30.42327 

-58.4060 

-40.5017 

79-47A 

46.06976 

-245.7439 

-196.6044 

11416U 

10.11341 

9.8689 

9.8849 

12553 

4.39037 

4.4278 

4.4316 

(a,  e,  i  are  the  same  as  those  in  Table  2-2.) 


Figure  2.2  SECULAR  CHANGES 
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2.7  Secular  Radial  Variation 

The  radial  distance  to  the  satellite  expressed  in  terms  of  orbital  elements  is 
r  =  a(l-ecosE)  (2.24) 

Expressing  cos  E  in  terms  of  a  series  expansion  of  derivatives  of  Bessel  coefficients 

cos  E  =  -e/2  +  2  £  -L  _£l_  J  (ne)  cos  nM  (2.25) 

n=i  n  d(ne)  n 

Substituting  equation  (2.25)  in  (2.24) 

2  00  i  j 

r  =  a+  -f-  -2ael  Fgg^Ja(ne)cosnM  (2.26) 

(see  Appendix  A  for  derivatives  of  Bessel  coefficients).  The  differentials  of  the 
independent  variables  a,  e,  and  M  are  defined  by 

da  =  Aa,  de  =  Ae,  dM  =  AM 


The  differential  of  the  dependent  variable  r  is  defined  by 

Ar  =  dr  =  tyfa  da  +  de  +  3r/dM  dM  (2.27) 

If  the  secular  change  in  radial  distance  is  required  and  only  the 
perturbation  due  to  the  gravitational  field  is  considered  then 


and 


A**sec  A®sec  ® 
Arscc  =  3f/dM  AMsec 


Using  (2.26) 


dr 


'dM  =  2aen?,  d(^)J°(ne)sinnM 


(2.28) 


From  either  (2.26)  or  (2.28)  it  is  seen  that  for  a  given  a  and  e,  there  is  no  secular 
change  in  r  for  a  given  M  +  2m7t  (m  is  an  integer).  Therefore, 


dr 


=  0 

M  +  2mrt 


Ar  ) 

sec  / 


=  o 


M  +  2m7t 


and 


From  orbit  to  orbit,  for  a  given  M,  there  is  no  secular  change  in  radial 
distance  if  the  only  force  acting  on  a  satellite  is  due  to  a  gravitational  field  with  axial 


symmetry. 

2.8  Freezing  Perigee  and  Eccentricity 

Low  altitude  near  circular  orbits  are  influenced  greatly  by  perturbations 

caused  by  the  Earth’s  oblateness.  The  magnitude  of  the  effects  of  the  perturbations 
depends  on  the  initial  values  of  certain  orbital  elements.  The  appropriate  selection  of 
semimajor  axis,  eccentricity,  inclination,  and  argument  of  perigee  can  minimize  the 
variations.  A  pertinent  example  consists  of  stopping  the  line  of  apsides  and  freezing 
the  perigee  point  [55-57].  If  the  periodic  change  in  eccentricity  is  considered  and 
equation  (2.18)  is  modified  to  include  the  periodic  term  of  J3  the  resulting  equations 
are: 


^dt  =  ~3/2  n^3  (R/p)3  0'e2)  sin  i  cos  co  (1  -  5/4  k)  (2.30) 

dW/dt  =  3nJ2  (%)2  (1  -  5U  k)  +  3n/2e  J3  (R/p)3  sin  03  sin  * 

x  [(1  -  5/4  k)  +  (35/4  cos2  i  -  cosec2  i)  e2  ]  (2.31) 

The  idea  of  this  type  of  frozen  orbit  is  to  initially  set  co  =  90°  causing  cle/clt 
due  to  J3  to  be  equal  to  zero.  Finding  the  root  to  equation  (2.31)  yields  the  value  of 
eccentricity  for  which  dw/dt  due  to  J2  and  J3  equals  zero.  These  values  of  e  and  co  will 
produce  a  perigee  location  which  remains  comparatively  constant.  The  changes  in  the 
average  argument  of  perigee  and  the  average  eccentricity  become  zero  when  co  =  90° 
and  e  approaches  some  determinable  small  value.  Cutting  (1978)  asserted  that 
subsequent  analysis  for  the  SEAS  AT- A  satellite  which  included  zonals  up  to  degree  21 
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lowered  the  appropriate  eccentricity  only  about  20%  [561. 

The  frozen  eccentricity  for  Landsat-5,  a  Sun  synchronous  orbit,  launched 
on  March  1,  1984,  is  approximately  0.0012.  The  averaged  argument  of  perigee  and 
eccentricity  will  then  oscillate  within  a  small  range  about  the  frozen  condition.  Fig.  2.3 
shows  the  evolutions  of  these  parameters  over  a  1 16  day  period  [55]. 

Some  other  satellite  orbits  which  would  satisfy  this  type  of  freezing  have 
the  following  elements 

a  (km)  i  (deg)  e 


7000  97.87  0.00103 

7188  98.57  0.00100 

7500  100.04  0.00094 


The  oscillations  of  to  and  e  can  be  observed  by  plotting  the  components  of 
the  eccentricity  vector;  the  result  is  a  circle  centered  on  the  frozen  condition  (see  figure 
2.4).  The  period  (T)  of  the  oscillation  is  related  to  the  secular  rate  of  the  change  in  to 
according  to 

T  _  2  7t 

A  toscc 

where  Aco^  is  obtained  from  equation  (2.18). 


i  i  r 


5  90.00  92.25  94.50  96.: 

T  OF  PERIGEE  (DEGREES) 

>F  e  AND  w  FOR  A  NEAR-FROZEI 
AT-5  ORBIT  [55] 


Figure  2.4  OSCILLATION  OF  e  AND  CD  ABOUT  FROZEN  POINT 
FOR  A  TYPICAL  NEAR-FROZEN  ORBIT  [55] 


CHAPTER  III 


VARIATIONS  DUE  TO  ATMOSPHERIC  DRAG 

3.1  Aerodynamic  Forces 

For  near-Earth  satellites,  excluding  the  Earth's  oblateness,  the 
perturbation  having  the  greatest  consequences  is  the  resisting  force  of  the  atmosphere. 
This  resisting  force,  called  drag,  is  especially  predominant  from  150  to  600  km.  When 
a  satellite  with  substantial  eccentricity  is  near  perigee  the  drag  tends  to  retard  its 
motion.  A  decrease  in  energy  causes  a  reduction  in  the  apogee  height;  this  is  called 
orbital  decay.  The  more  circular  the  orbit,  the  more  interminable  the  energy-draining 
effect  of  drag. 

The  main  effects  of  the  force  are  secular  decreases  in  the  orbital  elements 
a  and  e  which  in  turn  cause  the  orbit  to  contract.  Since  the  apogee  height  decreases 
while  the  perigee  height  remains  nearly  constant,  the  orbit  tends  to  become  more 
circular.  In  addition,  the  rotation  of  the  atmosphere  subjects  the  satellite  to  small  lateral 
forces  which  alter  slightly  the  orientation  of  the  orbital  plane.  The  result  is  a  small  but 
steadily  increasing  change  in  inclination  and  limited  periodic  changes  in  the  line  of 
apsides.  Small  changes  in  argument  of  perigee  are  caused  by  the  oblateness  of  the 
atmosphere  [4,  9,  58]. 

The  most  common  expression  for  the  drag  force  is 

D-^Co^m/P  Mv  (3-1) 

where  D  is  the  force  per  unit  of  mass,  m  is  the  mass  of  the  satellite,  CD  is  a 


dimensionless  drag  coefficient  which  is  empirically  determined,  and  A  is  the  effective- 
cross  sectional  area  of  the  satellite.  V  =  va  -  v,  where  v  is  the  satellite's  inertial  velocit} 
and  va  is  the  inertial  velocity  of  the  atmosphere,  /"represents  the  effect  of  atmospheric- 
rotation  on  drag  and  is  a  function  of  the  inclination  of  the  orbit.  Since  f  will  usually  be 
between  .9  and  1.1,  this  study  will  use  an  average  value  of  1.0.  For  details  on  the 
development  of  f  see  reference  58. 

The  ballistic  coefficient,  defined  by 


is  a  measure  of  the  ability  of  the  spacecraft  to  overcome  air  resistance.  For  typical 
spacecraft  /n  values  range  from  .01  x  10'6  to  .04  x  10'6  km2 /kg. 


Therefore, 


>!vi  v 


While  drag  is  acting  in  the  direction  opposite  to  the  satellite's  motion, 
there  are  forces  perpendicular  to  the  direction  of  motion.  These  forces  have 
components  passing  through  the  center  of  mass  which  cause  lift,  and  components 
which  produce  a  turning  moment  about  the  center  of  mass. 

The  turning  moment  is  normally  a  destabilizing  one  for  uncontrolled 
satellites.  Consequently,  the  resultant  value  of  the  variant  lift  force  will  be  zero.  For  a 
spinning  body,  there  are  moments  which  affect  the  orientation,  thereby  influencing  the 
lift  and  drag  forces.  The  precise  prediction  of  the  orbital  motion  of  elongated  bodies  of 
revolution  becomes  infeasible  in  the  absence  of  complete  information  about  orientation 
and  spin  everywhere  on  the  orbit.  For  these  reasons,  for  near-spherical  and  cylindrical 
satellites  with  a  length/diameter  ratio  greater  than  about  one,  the  forces  perpendicular  to 
the  direction  of  motion  can  be  ignored  [59]. 
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3.2  Density  Models 

The  feature  of  the  atmosphere  which  has  a  profound  effect  upon  the 
satellite  is  the  density  (p).  The  density  has  been  investigated  over  the  years  by  the 
study  of  the  rate  of  contraction  of  satellite  orbits  as  well  as  with  the  aid  of  instruments 
contained  in  the  satellites. 

Research  has  revealed  that  the  change  in  density  in  the  upper  atmosphere 
above  250  km  can  be  large  and  is  very  time-dependent.  Variation  in  density  occurs  as 
the  atmosphere  responds  to  solar  activity.  Fluctuations  in  density  include:  (1)  irregular 
diurnal  variations  resulting  from  ephemeral  solar  disturbances,  (2)  a  27-day  variation 
correlated  with  the  period  of  axial  rotation  of  the  Sun  with  respect  to  the  Earth,  (3)  an 
eleven  year  variation  corresponding  to  the  sunspot  cycle,  (4)  seasonal  variations 
caused  by  change  in  the  distance  from  the  Sun  and  the  inclination  of  the  Earth's  orbit 
to  the  Sun's  equator,  and  (5)  irregular  and  impermanent  variations  affiliated  with 
transient  geomagnetic  disturbances. 

Eliminating  the  time-dependent  variations  dictates  an  approximate  model 
of  the  atmosphere.  Important  assumptions  provide  analytical  representations  which 
allow  simplification  while  maintaining  reasonable  accuracy.  A  significant  simpli¬ 
fication  is  accomplished  by  assuming  that  the  atmosphere  is  spherically  symmetrical. 
Consequently,  the  atmospheric  density  is  a  function  of  the  distance  from  the  center  of 
the  Earth.  The  slight  inaccuracy  introduced  by  this  assumption  facilitates  the 
mathematical  formulation  necessary  for  determining  variations  in  orbital  elements. 

The  simplest  density  model  assumes  that  the  atmospheric  density 
decreases  exponentially  with  the  distance  from  the  center  of  the  Earth.  Assuming  that 
the  scale  height  (H)  is  constant  over  a  small  altitude  interval,  the  density  function  can 
be  written  as  [58,61  ]: 


p2  =  p,  expf(r,  -  r2)  /  H] 


(3.3) 


Consequently 


H  =  r‘  p2  (3.4) 

In  (Pi/P2) 

If  the  oblateness  of  the  atmosphere  is  to  be  considered,  the  same  density 
variation  (3.3)  can  be  used  as  a  function  of  the  altitude  above  the  Earth's  ellipsoidal 
surface.  The  altitude  is  determined  by  subtracting  the  radial  distance  to  the  surface  of 
the  Earth  from  the  radial  distance  of  the  satellite. The  radial  distance  to  the  surface  is 
found  using  [12] 

rE  =  R  [  1-  (f  +  3/2  f2  +...)  sin2  tp  +  3/2  f2  sin4  (p  -...]  (3.5) 

here 

f  =  flattening  ratio  of  the  Earth  =  V298.24 
<p  =  latitude 

R  =  equatorial  radius  =  6378.163  km 

Using  reference  or  cardinal  altitudes  (10  km  increments)  on  both  sides  of 
the  computed  altitude  (h  =  rE  -  r),  a  standard  atmospheric  table  such  as  the  U.S. 
Standard  Atmosphere,  1976  (see  Appendix  B)  can  be  used  to  compute  the  density  at 
the  two  radial  distances.  Exercising  equation  (3.4)  a  scale  height  can  be  computed  for 
that  particular  10  km  interval.  If  one  end  point  of  the  interval  is  now  used  as  position 
one  with  the  satellite  designated  as  position  two,  equation  (3.3)  can  be  employed  to 
determine  the  density  at  that  exact  altitude.  This  method  for  determining  density  is  best 
utilized  when  a  numerical  integration  scheme  is  used  to  continuously  update  the  values 
of  the  osculating  orbital  elements.  The  use  of  tabular  values  provides  the  accuracy  of  a 
dynamic  atmosphere. 

A  generalization  of  the  exponential  density  function  is  the  spherical 
power  density  function  model  which  assumes  the  scale  height  is  linear  [62] 


m 


(3.6) 
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where  t,  which  is  restricted  to  positive  integer  values,  and  s  are  disposable  parameters. 
For  coverage  of  all  reasonable  regions  averaged  over  an  1 1  year  cycle,  the  following 
values  result:  (1)  t  =  4,  (2)  rj  =  6498  km,  and  (3)  s  =  6456  km. 


3.3  Numerical  Computation  of  Element  Change 

To  examine  the  change  in  a  and  e  caused  by  drag  the  radial,  transverse, 
and  normal  force  components  (same  direction  as  fj,  f2,  and  f3  respectively  in  figure 
2.1)  are  expressed  in  terms  of  a  component  along  the  tangent  (fT)  and  a  component 
normal  to  the  orbit  (fN).  If  the  angle  which  the  tangent  to  the  orbit  makes  with  the 
transverse  component  is  y,  then 

fj  =  fT  sin  \j/  -  fN  cos  y 
f2  =  fT  cos  Y  -  fjsj  sin  Y 

These  values  for  fj  and  f2can  be  substituted  into  the  Lagrange  planetary 
equations  (2.5)  and  (2.6).  Since  the  force  normal  to  the  orbital  plane  is  extremely  small 
in  comparison  to  the  tangential  component,  it  can  be  ignored  for  this  approach.  This 
yields  a  modified  set  of  Lagrange  equations.  The  independent  variable  can  be  changed 
from  time  to  true  anomaly  using  two-body  relationships  and  subsequently  to  eccentric 
anomaly.  The  result  is 


da  -pa2(l+ e  cos  E)3^2 
^  P  (1-  e  cos  E)1/2 

de  __  -pa(l+  e  cos  E)1/2  (1-  e2)  cos  E 
^  P  (1-  e  cos  E)^2 

For  the  change  in  a  and  e  over  one  revolution 


(3.7) 

(3.8) 


Aa 


■  a_2J‘ 

P  o 


■  n  p(l+ecosE)3/2  dE 
>/2 


(3.9) 


Ae  =  -  J_[ 
r  j 


(1-  e  cos  E) 

2n  p  (1+  e  cos  E)'/2  (1- e2)  cos  E  dE 


PJo 


(1-  e  cos  E) 


h 


(3.10) 


Solv  ing  for  R  from  equation  (3.15)  and  substituting  back  into  (3.14) 

r  =  rp  (1-  f  sin2  ) 

(l-fsin2<Pp  ) 


(3.16) 


f 

I 

I 


I 

I 

I 


The  spheroid  is  now  defined  in  terms  of  latitude  of  the  initial  perigee  and 
latitude  of  the  satellite  at  any  position.  The  density  at  any  point  can  be  defined  in  terms 
of  the  density  at  the  perigee.  From  equation  (3.3) 


P2  =  Pp  exp  [  (  r '  r2  )  /  H  }  (3.17) 

Using  trigonometric  identities,  the  latitude  of  the  satellite  can  be 
expressed  as  a  function  of  the  position 

sin  tp  =  sin  i  sin  (to  +  v)  (3.18) 

Substituting  (3.16)  and  (3.18)  into  (3.17),  the  density  as  a  function  of 
radial  distance  and  true  anomaly  is 

p2  =  pP  exp  [(rp  -  r2)/H  +  c  cos  2  (to  +  v)  -  c  cos  2  co0  +  0(f2) 
where  c  =  */2  f  rp  sin2  i  /H  19) 

Expanding  exp(c  cos  2  (to  +  v))  as  an  exponential  series  expansion  and  transforming 
terms  of  true  anomaly  into  terms  containing  eccentric  anomaly  equation  (3.19)  can  be 
written  as  [58] 


p2  =  qexp{(a  e  cos  E  -  a)}[l+  c  cos  2  (to  +  E)  -  2ce  sin  2  (to  +  E) 

x  sin  E  +  c2/4(l+  cos4(to  +  E))  -  ce2/2(cos2to  +  2cos2(to  +  E) 
-  3  cos  2(to  +  E)  cos  2  E)  -  c2e  sin  4  (to  +  E)  sin  E 

+  O  (ce3,  c2e2)]  (3.20) 

where  q  =  pp  exp(rP/pj  -  c  cos  2too) 
too  =  initial  argument  of  perigee 


By  performing  the  following:  (1)  substituting  (3.20)  into  a  series 
expansion  of  the  integrands  (3.9)  and  (3.10),  (2)  utilizing  the  integral  representation  of 


the  Bessel  function  (i.e.  In  (z)  =  2/n  j02n  exp  (z  cos  A)  cos  n  A  dA,  where  A  is  an 
angle),  and  (3)  eliminating  periodic  terms  the  following  represent  changes  over  one 
revolution 

Aa^  =  -27ia2/p  q  expCa/H)[I0+  2e  I,+  3/4  e2  (I0+  I2)  +  e3/4  (3  I1+  I3) 
+  0(e4,ce3)]  (3.21) 

Aesec  =  -2ira2/p  q  exP  Ca/H)[V  %  (Io+  I2>  +  e2/8  (-5  Ij+  I3) 
-e3/16(5I0  +  4I2-I4)  +  O(e4)]  (3.22) 

Substituting  equations  (3.21)  and  (3.22)  into  (3.11) 

=  -2tta2/p  q  exp('a/H)  [IQ- 1,-  e/2  (3  IG  -  4  Ij  +  I2) 

+  e2/g  (6  10-  1 1  I,+  6 12  - 13)+  e3/16  (-7  I0+12  Ir  8  I2+  4 13  - 14) 
+  0(e4)]  (3.23) 


The  argument  of  the  Bessel  functions  is  z  =  ae/j^.  I0  and  Ij  can  be  estimated  using  an 
asymptotic  expansion  or  polynomial  approximation  while  I2, 13,  and  I4  can  be  found 
using  the  recurrence  relation 

In,l(z)  =  In-l<Z)'2n/zIn(Z) 

For  the  complete  development  of  the  change  in  the  other  orbital  elements 
over  one  revolution  see  reference  58. 

Ai  =^2w(-^j/2  exP  ("V  sin  j  H0  +  (1  +  c)I2  cos  2w 

-4el  i  cos2 co  +  y  c  (Iq  +  I4  cos  4co)  +  0(c2  e2  )]  (3.24) 

where  w  =  angular  velocity  of  the  atmosphere 

and  /  =  (1  -  cosi)2 

vP 

Substituting  the  half-angle  formula  for  cos2co  and  eliminating  the  periodic  terms 

Ai«c.^“'(-f^)/2exp(-a/H)sini[I0-2cIl)  (3.25)  jj 


If  (3.25)  is  written  as 

Aisec  =  -S  [1  o  -  2el]  ] 

it  is  evident  that  for  all  values  of  a,  e,  and  i,  S  >  0.  Therefore,  the  inclination  angle 
decreases  as  long  as  I0  >  2el!.  Using  a  polynomial  approximation  [63] 

I0(z)  =  1  +  .24999985  z2  +  .01562519  z4  +  ... 

Ij(z)  =  z[V2  +  .062499978  z2  +  .00260419  z4  +  ...] 
it  is  clear  that  for  all  feasible  values  of  eccentricity  for  near  Earth  satellites  ID  >  2el,. 

Examining  a  worst  possible  situation  demonstrates  that  the  relationship  is 

satisfied 


a  =  7000  km 
H  -  72.6 


I0  =  26.868438 


2el.  -  3.8241799 


e  =  .06 


The  other  two  orbital  elements,  Q  and  to,  experience  no  secular  change 


due  to  drag. 


CHAPTER  IV 


i 

MINIMUM  ALTITUDE  VARIATION  ARCS  | 


4.1  Requirement 

A  satellite  which  experiences  minimum  altitude  variation  over  a 
circumscribed  latitude  range  can  be  very  useful  for  scientific  research  and  application 
experiments.  Motivation  for  achieving  and  maintaining  this  category  of  frozen  orbit 
includes:  (1)  study  of  the  flattening  and  other  dynamic  properties  of  the  atmosphere, 
(2)  photo  resolution  optimization  of  preset  optical  systems,  and  (3)  minimization  of 
necessary  geometric  corrections  to  satellite-produced  images. 

A  good  example  of  the  necessity  of  this  type  of  arc  is  the  SEASAT-A 
satellite,  launched  on  June  26,  1978,  to  provide  a  global  ocean  dynamics  monitoring 
system.  The  instrument  used  to  obtain  radar  imagery  of  ocean  waves,  ice  and  snow 
cover,  and  land  surfaces  is  a  Synthetic  Aperture  Radar  (SAR)  which  has  a  resolution 
of  25  meters.  Constraints  dictate  minimum  altitude  variation  over  SAR  stations  which 
are  all  located  in  the  Northern  Hemisphere  [56]. 

To  date,  very  little  literature  exists  on  the  subject  of  minimum  altitude 
variation  arcs.  Some  work  was  accomplished  by  Kalil  [64]  as  he  studied  the 
first-order  effects  of  oblateness  while  neglecting  entirely  atmospheric  drag.  Since  it  is 
near-Earth  satellites  for  which  the  requirement  for  minimum  altitude  deviation  may 
exist,  it  is  essential  to  include  the  effects  of  atmospheric  drag  to  realize  a 
comprehensive  study. 
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Altitude  deviations  are  best  minimized  by  essentially  modeling  the  shape 
of  the  Earth  over  the  latitude  range  of  interest.  The  proper  combination  of  semimajor 
axis,  eccentricity,  inclination,  and  argument  of  perigee  produces  the  desired  arc. 


4.2  Latitude  Coverage  and  Eccentricity  Required 

To  determine  the  minimum  altitude  variation  several  approaches  can  be 
exercised.  One  such  approach  involves  the  continuous  computation  of  the  altitude 
along  the  orbit  at  intervals  defined  by  a  predetermined  step  size.  The  difference 
between  the  continuously  updated  altitude  and  the  initial  altitude  is  computed  until  the 
maximum  allowable  limit  is  reached.  This  procedure  becomes  a  tradeoff  between 
eccentricity  and  the  latitude  coverage  attainable. 

Referring  to  figure  4.1:  rs  is  the  radial  distance  from  the  center  of  the 
Earth  to  the  satellite,  rE  is  the  distance  from  the  center  to  the  surface  of  the  Earth,  and  h 
is  the  altitude  above  the  surface  of  the  Earth.  The  starting  position  is  designated  as  1 
while  2  represents  the  position  at  an  interval  or  multiple  of  intervals  of  time  later.  Each 
time  an  interval  is  covered  the  new  state  becomes  position  2. 


Figure  4.1  RADIAL  DISTANCES  AND  ALTITUDES 


The  altitudes  corresponding  to  the  respective  positions  are 


hl  - h,  -  % 

(4.1) 

h2  =  rs2  “  rE2 

(4.2) 

Subtracting  the  above  equations 

h,  -  h2  =  ('»,  -  h2)  ■  (%  -  <e2) 

(4.3) 

Substituting  equation  (2.26)  for  rs  and  equation  (3.5)  for  rE,  assuming  a 
secularly  constant  a  and  e  (no  drag),  and  using  (4.3)  to  include  0(e3)  for  e  <  .1 
h,  -  h2  =  a  [(cos  M2  -  cos  Mj,e  +  (cos  2M2  -  cos  2Mj)e2 
+  3/8  (3  cos  3M2  -  3  cos  3Mj  -  cos  M2  +  cos  Mj)e3] 

-  R  [(f  +  3/2  f^Xsin2^  -  sin2<pj)  +  3/2  ^(sin4^,  -  sin4cp2)] 

(4.4) 

There  are  two  ways  in  which  equation  (4.4)  can  be  employed.  The  first 
approach  is  to  compute  the  latitude  range  that  can  be  achieved  prior  to  exceeding  the 
maximum  allowable  deviation  from  the  initial  altitude  given  an  a  and  e.  The  second 
way  is  to  start  with  the  semimajor  axis  and  latitude  range  desired  and  then  determine  a 
compatible  e  which  permits  |  hj  -  h2  |  to  be  a  minimum.  This  method  requires  an 
iterative  scheme  where  I  h j  -  h2 !  =  hj  2  is  defined  as  the  performance  index  [65].  For  a 
minimum  to  exist  over  the  range  of  possible  values  of  eccentricity, 

Ahj_2  (e)  =  hj_2  (e  +  Ae)  -  h12  (e)  >  0 
where  Ae  is  a  small  incremental  change. 

The  latitude  can  be  computed  at  the  end  of  each  interval  in  the  following 

manner: 

(1)  An  initial  satellite  radial  distance  and  mean  anomaly  is  established. 

(2)  Recall  Aasec  =  Aesec  =  0  in  the  absence  of  drag.  Using  equation  (2.28) 
and  the  secular  equation  (2.19)  for  the  rate  of  change  of  mean  anomaly 
per  second,  the  radial  distance  can  be  computed  at  each  new  position. 


dr  =  Ar  =  ^r/dM  dM 
rs2  =  rSj  +  Ar 

(3)  The  updated  mean  anomaly  is 

M2  =  Mj  +  AM 

(4)  From  Kepler's  equation,  M  =  E  -  e  sin  E,  a  method  of  successive 
approximations  can  be  employed  to  determine  eccentric  anomaly. 

(5)  True  anomaly  for  the  elliptic  orbit  can  be  calculated  from  the  relation 

(6)  Given  an  argument  of  perigee  co,  the  argument  of  latitude  is 

u  =  v  +  to  (4.5) 

(7)  From  spherical  trigonometry 

sin  cp  =  sin  i  sin  u  (4.6) 

(8)  All  variables  are  now  known  for  substitution  into  equation  (4.4)  to 
compute  the  deviation  in  altitude  from  the  initial  position  to  the  next 
position  of  the  satellite. 

(9)  If  e  is  the  predetermined  maximum  allowable  altitude  deviation,  the 
process  is  repeated  until  |  hj  -  h2 1  >  e 


The  radial  distance  measured  from  the  center  of  the  Earth  to  the  satellite 


in  terms  of  true  anomaly  is 


a (1-  e2) 
1+  e  cos v 


Equation  (3.5)  to  O(f)  is 


rE  =  R  [1-  f  sin2  cp] 


(4.9) 


substituting  (4.5)  and  (4.6)  in  (4.8) 
rE  =  R  [1-f  sin2  i  sin2  (v  +  to)] 

The  altitude  above  the  surface  of  the  Earth  is 

h  =  r-rE  (4.10) 

Substituting  (4.7)  and  (4.9)  in  (4.10),  assuming  no  drag  (i.e.  a  and  e 
are  secularly  constant),  taking  the  derivative  of  altitude  with  respect  to  true  anomaly, 
and  simplifying  with  trigonometric  identities  the  result  is 

^  =  “111^15111^+ Rf  sin2  i  sin  2(  v  +  co)  (4.11) 

av  (1+  e  cos v)2 

From  (4.1 1)  it  is  obvious  that  the  change  in  altitude  is  a  function  of  the 
shape  of  the  Earth,  the  size  and  shape  of  the  orbit,  and  the  orientation  of  the  orbital 
plane.  For  a  given  a  and  e  the  rate  of  change  is  highly  dependent  upon  the  argument  of 
perigee  and  inclination  angle.  For  a  given  co  the  change  in  altitude  with  respect  to  the 
change  in  true  anomaly  is  a  periodic  function  with  period  2k. 

The  next  series  of  graphs  demonstrate  the  effect  of  argument  of  perigee, 
eccentricity,  semimajor  axis,  and  inclination  on  the  rate  of  change  of  altitude  as 
determined  by  equation  (4.11).  A  comparison  of  figures  4.2  and  4.3  shows  that  a 
change  of  500  km  in  semimajor  axis  has  little  effect  on  the  rate  of  change  of  altitude  for 
a  given  eccentricity,  inclination,  and  argument  of  perigee.  The  maximum  change  in  rate 
is  ±  5  km/rad. 

A  change  in  eccentricity  from  .01  to  .0024  causes  a  large  change  for  a 
given  semimajor  axis,  inclination,  and  argument  of  perigee.  From  figures  4.3  and  4.4 
the  change  in  the  rate  is  seen  to  be  as  high  as  50  km/rad.  Figures  4.4,  4.5,  and  4.6 
show  that  for  a  constant  a,  e,  and  co,  as  the  inclination  angle  increases  to  90°  the  rate  of 
change  of  altitude  increases.  From  i  =  45°  to  i  =  90°  the  change  in  the  rate  is  about 


lOkm/rad.  Study  of  the  graphs  reveals  that  the  minimum  rate  of  altitude  change 
(d^/dv  =  0)  occurs  at  or  around  perigee/apogee.  Regardless  of  the  argument  of  perigee 
as  the  eccentricity  increases  the  minimum  rate  progresses  closer  to  the  perigee/apogee. 
Therefore,  the  minimum  altitude  variation  occurs  more  closely  to  perigee/apogee. 

Zero  rate  of  change  always  occurs  at  the  perigee  when  it  is  located  on  the 
equator  (i.e.  co  =  0°).  As  the  eccentricity  decreases  the  argument  of  perigee  causes  the 
value  of  true  anomaly  for  which  ^/dv  =  0  to  shift  away  from  perigee/apogee.  Figure 
4.5  shows  that  this  displacement  is  as  much  as  70°. 

The  latitude  range  over  which  the  minimum  altitude  deviation  occurs  can 
be  determined  easily  from  the  figures.  All  three  values  of  i,  v,  and  to  corresponding  to 
the  minimum  rate  can  be  obtained  from  the  graph  and  equation  (4.6)  used  to  determine 
the  respective  latitude. 

Data  contained  in  figures  4.7,  4.8,  and  4.9  include  oblateness  effects. 
Figure  4.7  demonstrates  that  the  minimum  altitude  a  satellite  obtains  in  an  orbit  does 
not  necessarily  occur  at  the  same  latitude  at  which  perigee  occurs.  The  latitude  of  the 
perigee  tpp  on  this  plot  is  50.77°N.  Depending  on  the  eccentricity  the  minimum  altitude 
deviates  from  perigee.  As  the  eccentricity  increases  the  secular  change  in  the  altitude 
over  the  latitude  range  of  interest  increases. 

Figures  4.8  and  4.9  show  a  magnified  view  of  the  effect  of  argument  of 
perigee  over  a  latitude  range  30°  to  60°N.  The  latitude  range  of  interest  can  vary 
depending  upon  the  objective  of  the  satellite  mission.  Deviation  is  measured  against  the 
altitude  at  the  start  of  the  arc.  In  this  case,  study  of  the  entire  latitude  arc  of  30°  shows 
that  if  the  orbit  has  an  eccentricity  of  .005  and  initial  argument  of  perigee  of  66°  the 
change  in  altitude  is  +  6km.  If  instead  to  =  80°  the  deviation  over  30°  of  latitude  is  less 
than  1  km.  By  proper  adjustment  of  argument  of  perigee  the  altitude  deviation  over  any 
latitude  range  can  be  minimized  for  a  given  semimajor  axis  and  eccentricity. 
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Starting  with  the  combination  of  (4.5)  and  (4.6) 
sin  <p  =  sin  i  sin  (v  +  co) 

Using  the  trigonometric  identity  for  the  sum  of  angles 
sin  <p  =  sin  i  [sin  v  cos  co  +  cos  v  sin  co] 

Since  v  =  0  at  the  perigee  latitude  <pp 


sin  cpp  =  sin  i  sin  co  (4. 12) 

The  derivative  of  (4.12)  with  respect  to  time  is 

cos  <pp  <pp  =  sin  i  cos  co  cb  +  sin  co  cos  i  tydt 


Recalling  <^Vdt)sec  =  0 

The  secular  change  in  the  perigee  latitude  is 


<pn  =  sin  i  cosco  co 
p  cos  <Pp 


(4.13) 


From  (4. 1 3)  it  is  obvious  that  there  are  two  ways  to  freeze  the  perigee 
latitude.  The  first  is  by  achieving  an  argument  of  perigee  of  90°.  As  discussed  in 
section  2.8  when  co  =  90°  the  de/dt  due  to  the  J3  term  is  zero.  The  eccentricity  can  then 
be  computed  so  that  d°tydt  due  to  the  J2  and  J3  terms  is  zero.  The  second  way  is  to 
establish  an  orbit  having  an  inclination  angle  that  causes  the  change  in  co  to  be  zero  for 
a  given  eccentricity  and  semimajor  axis.  As  discussed  in  section  2.5  to  0(J2)  this 
occurs  when  i  =  63.44°. 
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Figure  4. 10  THE  LATITUDE  ARC 


If  the  latitude  arc  is  divided  into  N  equal  segments  fixed  by  a 
predetermined  time  increment  the  altitude  change  at  the  end  of  each  interval  can  be 


determined  in  the  following  manner: 

hi  =  VrEl 

h2  =  rs,  +  ^1,2  -  rE, 


(4.14) 

(4.15) 


Subtracdng(4.15)  from  (4.14) 

hj  -  h2  =  Ah12  =  rEz  -  rEj  -  Ar12 

For  the  next  interval 

Ah2,3  =  rE3  '  rE2  '  Ar2,3 
In  general  terms 

AhN-l,N  =  rEN'rEN.1  -  ^N-l.N 


Using  (2.27)  for  ArN1  N  the  altitude  difference  over  each  interval  on  the 


arc  is  given  by 
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^"Ahl,2  ^ 

f%  -  % 

/%Aa  1.2  +  +|JiAM1.2^ 

Ah2>3 

= 

%  ;  % 

- 

^Af2.3  +  ^Ae2,3  +|^A;M2.3 

ls.AhN-l,N^ 

Ven  -  >en-J 

^  ^AaN-l,N+  ^§AeN-l,N  +  |ftAMN-l,^ 

(4.16) 

The  value  of  the  orbital  elements  can  be  computed  using  an  Euler  type 
scheme.  From  the  initial  values  of  a,  e,  M,  and  E  at  the  beginning  latitude  the  change 
in  mean  anomaly  over  the  first  interval  can  be  determined  from  (2.19).  The  mean 
anomaly  at  the  start  of  the  next  interval  is 
M2  =  +  AMj  2 

E2  can  be  calculated  from  M2  using  an  iterative  process.  Employing  (3.9)  and  (3.10) 
the  change  in  semimajor  axis  (Aaj  2)  and  eccentricity  (Ae1>2)  due  to  drag  can  be 
computed  using  Ej  and  E2  as  the  limits  of  integration.  There  is  no  secular  change  in 
these  two  elements  due  to  the  gravitational  field.  However,  secular  changes  due  to 
gravity  occur  over  the  interval  to  Q  and  0)  and  can  be  calculated  if  desired.  Therefore, 

&2  =  aj  +  Aaj  2 

e2  =  ei  +  Ae1(2 

Q2  =  +  (dfl/dt)  Atj  >2 

G>2  =  CDj  +  (d^/dt)  Atj  2 

With  the  new  values  of  a  and  e  an  updated  value  of  mean  anomaly 
change  AM2  3  can  be  determined.  Then 

M3  =  h42  +  AM2  3 

The  process  is  repeated  until  N  intervals  have  been  completed  and  the 


exit  latitude  is  reached. 


Using  steps  5  and  6  on  page  50  and  equation  (4.9)  rE  can  be  calculated  at 
the  end  of  each  interval.  All  information  is  now  available  for  substitution  into  (4.16). 
The  equation  can  be  to  compute  the  altitude  deviation  over  each  interval  as  well  as 
the  maximum  positive  and  negative  deviation  over  the  entire  arc.  Additionally,  the 

minimum  amount  of  total  altitude  deviation  over  the  arc  occurs  when 

N 

^  (AIin-i,n) 

n=z 

is  minimized.  In  all  cases  it  becomes  a  matter  of  selecting  the  appropriate  a,  e,  and  to  to 
minimize  the  deviation. 

For  this  study,  a  fourth-order  Runga-Kutta  method  was  used  in 
combination  with  a  tabular  density  model  (see  Appendix  B)  for  the  integration  of  (3.9) 
and  (3.10).  Three  different  step  sizes  of  3,  5,  and  10  seconds  were  exercised;  the 
results  varied  very  little.  A  fraction  of  the  findings  are  presented  in  Tables  4-1  through 
4-6. 

Max  +  is  the  maximum  positive  altitude  deviation  for  the  latitude  range 
studied  while  Max  -  is  the  maximum  negative  altitude  deviation.  Both  deviations  are 
measured  against  the  altitude  at  the  initial  latitude,  hj  -hN  is  the  altitude  difference 
between  the  initial  altitude  and  the  altitude  N  intervals  later  at  the  exit  latitude.  In  Tables 
4-5  and  4-6  Aa  is  the  change  in  semimajor  axis  over  the  latitude  range  while  Ae  is  the 
change  in  eccentricity. 

Table  4-2  shows  a  comparison  of  minimum  altitude  deviation  for  a  very 
small  eccentricity  (.0022)  versus  a  larger  eccentricity  (.01).  Minimum  deviation  for  the 
higher  eccentricity  is  approximately  2.7  km  versus  .36  km  for  the  lower  e. 

Tables  4-1,  4-2,  and  4-3  demonstrate  the  effect  of  eccentricity  and 
argument  of  perigee  on  minimum  altitude  deviation  for  a  given  semimajor  axis  and 
inclination.  Tables  4-4  and  4-5  show  the  effect  of  semimajor  axis  on  altitude  deviation 


over  a  latitude  range  given  an  eccentricity,  inclination,  and  argument  of  perigee.  The 
effect  of  drag  is  to  decrease  the  maximum  positive  deviation  while  increasing  the 
maximum  negative  deviation.  Table  4-6  shows  a  comparison  of  data  determined  from 
the  same  orbital  elements  but  different  ballistic  coefficients  and  different  integration 
step  sizes.  As  (3  increases  the  maximum  positive  altitude  deviation  decreases  while  the 
maximum  negative  deviation  increases. 

The  data  in  Table  4-3  correlates  favorably  with  figure  4.5.  The  minimum 
altitude  deviation  for  e  =  .0024,  a  =  7500  km,  and  i  =  63.44°  occurs  when  to  =  69°. 
Over  the  latitude  range  of  30°-60°N  the  deviation  is  approximately  .18  km.  A  review 
of  figure  4.5  indicates  that  when  to  =  69°  the  minimum  rate  of  altitude  change  is 
experienced  between  35°  and  77°  true  anomaly  which  corresponds  to  latitudes  60°N 
and  30°N  respectively.  Examination  shows  that  over  this  range  the  ^/dv  curve  rides 
the  zero  value  for  approximately  15°. 

The  data  contained  in  the  tables  that  follow  include  secular  effects  only. 
The  variational  equations  of  the  averaged  elements  were  presented  without  having 
included  the  short-periodic  variations.  The  effects  contributed  by  the  short-periodic 
variations  can  be  studied  by  exercising  the  corresponding  equations  developed  in 
reference  54.  These  equations  can  be  used  with  (4. 16)  to  predict  the  altitude  deviation 
due  to  periodic  variations.  The  procedure  outlined  in  section  4.5  would  then  be 
accomplished  to  determine  the  combination  of  orbital  elements  that  would  produce 
minimum  altitude  deviation  over  the  desired  latitude  range. 


Table  4-2  ALTITUDE  DEVIATION  OVER  LATITUDE  RANGE  30°  -  60°N 

e  =  .0022  and  .01 
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a  =  7000  km 


i  =  63.44° 

e 

co(deg) 

Max  + 

(km) 

Max  - 

(km) 

hrhN 

(km) 

.0022 

50.0 

.4790 

.1356 

.1356 

.0022 

53.5 

.3634 

.3283 

.3283 

*  .0022 

54.0 

.3555 

.3527 

.3527 

.0022 

54.5 

.3283 

.3941 

.3941 

.0100 

62.0 

4.4216 

1.9143 

-4.4216 

.0100 

64.0 

2.7079 

2.5605 

-2.7079 

*  .0100 

64.5 

2.2976 

2.6576 

-2.2976 

.0100 

65.0 

1.9001 

2.8653 

-1.9001 

a  = 

i  = 

Table  4-3 

7500  km 

63.44° 

e 

ALTITUDE  DEVIATION  OVER  LATITUDE  RANGE  30°  -  60°N 

.002  <e<. 0032 

Max  +  Max  -  hj-hN 

to(deg)  (km)  (km)  (km) 

* 

.0020 

49.5 

.4150 

.3798 

.3798 

* 

.0023 

65.5 

.2197 

.1736 

.1322 

* 

.0024 

69.0 

.1632 

.1781 

.0944 

* 

.0025 

71.5 

.1395 

.1298 

-.0150 

* 

.0026 

74.0 

.1055 

.1323 

-.0771 

* 

.0027 

76.0 

.1796 

.1096 

-.1796 

* 

.0028 

78.5 

.1476 

.1956 

-.1476 

* 

.0032 

85.0 

.3059 

.3332 

-.3059 

Table  4-4  ALTITUDE  DEVIATION  OVER  LATITUDE  RANGE  25°-58°N 

No  Drag 


i  =  63.44° 

to  =  60.0° 

e  =  .0024 

a 

Max  + 

Max  - 

hrhN 

(km) 

(km) 

(km) 

(km) 

7000 

.1136 

6950 

.0987 

6900 

.0854 

6850 

.0819 

6800 

.0710 

6750 

.0614 

6700 

.0530 

6650 

.0391 

6600 

.0332 

6550 

.0281 

.1770 

.0111 

.2391 

.0752 

.3019 

.1491 

.3570 

.2269 

.4216 

.3021 

.4869 

.3884 

.5529 

.4644 

.6262 

.5474 

.6934 

.6329 

.7615 

.7097 

Table  4-5  ALTITUDE  DEVIATION  OVER  LATITUDE  RANGE 

With  Drag 


i  =  63.44° 

co  =  60.0° 

e  =  .0024 

(3  =  .04  x  10^  km2/kg 

a 

Aa 

Max  + 

Max  - 

(km) 

(km) 

Ae 

(km) 

(km) 

7000 

.0001 

.0000000 

.1136 

.1772 

6950 

.0002 

.0000000 

.0986 

.2395 

6900 

.0005 

.0000000 

.0853 

.3028 

6850 

.0010 

.0000000 

.0816 

.3591 

6800 

.0024 

.0000001 

.0703 

.4264 

6750 

.0056 

.0000003 

.0598 

.4987 

6700 

.0144 

.0000007 

.0493 

.5845 

6650 

.0413 

.0000020 

.0303 

.7223 

6600 

.1469 

.0000071 

.0090 

1.0760 

6550 

1.0295 

.0000486 

.0000 

3.8963 

25°-58°N 


hrhN 

(km) 


.0108 

.0757 

.1502 

.2295 

.3080 

.4028 

.5021 

.6599 

1.0589 

3.8963 


Table  4-6  ALTITUDE  DEVIATION  OVER  LATITUDE  RANGE  30°-60°N 

With  Drag 
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i  =  63.44° 

a 

(km) 

co  =  76.0° 

Aa 

(km) 

e  =.003 

Ae 

(3  =.01xl0'6 
Max  + 

(km) 

Interval  Size 

Max  - 

(km) 

=  5  sec. 

hrhN 

(km) 

6800 

.0007 

.0000001 

.2526 

.0646 

-.2526 

6750 

.0016 

.0000001 

.1632 

.1179 

-.1632 

6700 

.0041 

.0000004 

.0950 

.1671 

-.0950 

6650 

.0118 

.0000011 

.0647 

.2324 

.0089 

6600 

.0429 

.0000041 

.0501 

.3342 

.1470 

6550 

.2523 

.0000245 

.0156 

.7736 

.7189 

P  =.04xl0"6 

Interval  Size 

=  5  sec. 

6800 

.0026 

.0000002 

.2491 

.0667 

-.2491 

6750 

.0063 

.0000006 

.1624 

.1154 

-.1624 

6700 

.0164 

.0000016 

.0745 

.1817 

-.0721 

6650 

.0476 

.0000046 

.0548 

.2790 

.0779 

6600 

.1759 

.0000170 

.0314 

.5335 

.4180 

6550 

1.2844 

.0001244 

.0000 

2.8811 

2.8811 

P=.01xI0-6 

Interval  Size 

=  10  sec. 

6800 

.0007 

.0000001 

.2448 

.0724 

-.2448 

6750 

.0016 

.0000001 

.1462 

.1179 

-.1462 

6700 

.0041 

.0000004 

.0783 

.1670 

-.0783 

6650 

.0117 

.0000011 

.0579 

.2392 

.0158 

6600 

.0429 

.0000041 

.0501 

.3341 

.1469 

6550 

.2481 

.0000241 

.0119 

.7765 

.7298 

To  maintain  the  orbit  and  freeze  the  apogee  and  perigee  height 
maintenance  maneuvers  must  be  performed.  One  way  of  performing  the  maneuvers  is 
to  allow  the  orbit  to  decay  until  one  of  the  elements  has  changed  an  amount  equal  to  or 
greater  than  a  prescribed  tolerance  for  that  element  Another  way  is  to  apply  corrective 
propulsion  at  perigee  and  apogee  each  orbit.  In  this  study  the  second  approach  is  used 
to  provide  the  propulsion  for  counteracting  drag  effects.  In  both  cases  a  two-impulse 
maneuver  is  initiated  which  transfers  the  satellite  back  to  the  original  orbit.  The  transfer 
is  virtually  coplanar  since  the  atmospheric  velocity  is  small  compared  to  the  satellite 
velocity  and  consequently  the  majority  of  the  perturbing  force  is  in  the  plane  of 
motion. 

To  conserve  propellant  a  minimum  total  AV  is  of  particular  interest.  The 
associated  maneuver  which  corresponds  to  the  least  amount  of  energy  consumed  is  the 
Hohmann  transfer.  This  type  of  transfer  represents  the  minimum  total  AV  maneuver 
for  cases  where  the  ratio  of  large  to  small  orbit  radius  is  less  than  1 1 .8  [5 1 ,66]. 

The  AV  required  to  maintain  the  apogee  height  and  then  the  perigee  must 
be  computed  separately  and  then  added  together  to  determine  the  total  impulse 
required.  The  amount  of  AV  that  is  depleted  as  the  satellite  moves  from  the  perigee  to 
apogee  should  be  added  at  the  perigee.  The  amount  of  AV  that  is  dissipated  going  from 
apogee  back  to  the  perigee  should  be  inserted  at  the  apogee. 

Using  the  vis-viva  integral  the  total  AV  required  can  be  computed  in  the 


following  manner 


perigee 


( 1 )  no  drag  case 

(2)  actual,  if  no  corrective 
propulsion  (o  <  v  <  180) 


Figure  4. 1 1  APOGEE  CHANGE 
DUE  TO  DRAG 


The  velocity  at  apogee  if  drag  were  not  present 


Since  drag  affects  the  radius  and  apogee  height 


__M _ 

ar  AaA 


Therefore 


IV1A-V2A  ^  perigee 


(3)  =  actual,  if  no  corrective 
propulsion  (180  <v  <  360°) 


Figure  4.12  PERIGEE  CHANGE 
DUE  TO  DRAG 


apogee 


68 


The  velocity  at  perigee  if  drag  were  not  present 


Since  drag  affects  the  radius  and  perigee  height 


2p 


VArP 


J± _ 

3j+  Aap 


Therefore 


AV 

@ apogee 


And 


AV  =  AV  .  + 

Total  (gjpengee 


AV 

@ apogee 


Table  4-7  shows  the  AV  required  to  maintain  the  orbit  given  an  initial  set 
of  orbital  elements.  Al,  El,  RIP,  and  R1 A  are  the  semimajor  axis,  eccentricity,  radius 
of  perigee,  and  radius  of  apogee  for  an  unperturbed  orbit.  A2,  E2  and  R2A  represent 
the  orbital  parameters  that  exist  at  apogee  if  drag  is  not  compensated  for  at  the  perigee. 
A3,  E3,  and  R3P  represent  the  orbital  parameters  that  exist  at  perigee  if  drag  is  not 
compensated  for  at  the  apogee.  As  would  be  expected  the  lower  the  semimajor  axis  the 
more  the  AV  required  to  maintain  the  orbit. 

As  is  seen  from  Tables  4-4  and  4-5  the  smaller  the  semimajor  axis  the 
more  the  drag  affects  the  altitude  deviation  over  the  latitude  range.  Corrective 
propulsion  becomes  essential  to  maintaining  the  minimum  altitude  variation  arc. 


6525  6508.12  6496.84  .00311  .00085  6501.51  6548.49  6528.41  6491.29 


CHAPTER  V 


SUN  SYNCHRONOUS  ORBITS 

5.1  Sun  Synchronous  Orbits 

Since  the  early  1960s  the  Sun  synchronous  orbit  has  been  used  for  most 
every  meteorological  satellite  as  well  as  Earth  Resources  Technology  satellites.  Data 
from  satellites  in  this  type  of  orbit  will  generate  most  of  the  Earth  resources  data  in  the 
1985*1995  time  period  [67].  This  familiar  orbit  uses  the  precession  of  the  orbit  line-of- 
nodes  caused  by  the  Earth's  oblateness  to  maintain  a  fixed  angular  orientation  of  the 
orbit  plane  relative  to  the  Sun.  The  required  precession,  which  depends  upon  the 
inclination  and  semi-latus  rectum,  is  approximately  .9856°/day.  The  Sun  synchronous 
orbit  condition  occurs  when  the  orbit  nodal  motion  is  equal  to  the  motion  of  the  Earth 
about  the  Sun.  The  orbital  orientation  is  effectively  frozen  with  respect  to  the 
Earth-Sun  line.  By  careful  execution  of  final  insertion  into  the  desired  orbit  the  angle 
between  the  orbit  node  and  the  Sun's  right  ascension  will  remain  nearly  constant. 

The  Sun  synchronous  orbit  is  particularly  suitable  for  Earth  observation  and 
remote  sensing  applications  for  the  following  reasons:  (1)  the  satellite  traverses  each 
latitude  at  the  same  local  time;  the  consequence  is  similar  ground  lighting  conditions  on 
each  pass,  (2)  the  majority  of  the  Earth's  surface  can  be  mapped  with  near  north-south 
contiguous  swaths  in  a  fixed  period  with  repeatability,  and  (3)  the  average  spacecraft 
solar  array  angle  of  incidence  to  the  Sun  remains  within  defined  limits  ensuring 
continuous  electric  power. 
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5.2  Sun  Synchronous  Degradation 

Any  perturbation  which  alters  the  orbit  altitude  or  the  orbit  inclination 
promotes  the  loss  of  Sun  synchronism.  One  primary  source  of  synchronism 
degradation  is  atmospheric  drag.  The  principal  effect  of  drag  is  to  reduce  altitude  and 
thereby  cause  period  reduction.  This  orbit  decay  separates  the  corresponding  tracks  of 
two  successive  repeat  cycles.  Consequently,  the  tracks  are  no  longer  exactly  overlaid. 
The  study  presented  in  reference  69  reports  that  for  a  satellite  with  an  average  ballistic 
coefficient  initially  at  678  km  the  loss  in  altitude  after  three  months  is  .19  km.  The  loss 
in  altitude  induces  a  nodal  sideslip.  The  sideslip  results  in  a  loss  of  in-orbit  phasing  of 
approximately  2.6  minutes  which  corresponds  to  an  eastward  movement  in  nodal 
crossing  of  72.267  km.  If  corrections  are  made  every  three  months,  Sun  synchronism 
may  be  considered  as  varying  ±1.3  minutes  about  a  mean  node;  an  insignificant 
disturbance.  The  adjustment  can  be  achieved  by  the  two-impulse  Hohmann  transfer 
orbit.  The  AV  required  every  three  months  would  be  about  .0914  m/sec. 

Another  source  of  degradation  is  the  gravitational  force.  A  1975  Goddard 
Space  Flight  Center  study  indicated  that  even  zonal  harmonics  through  order  4  ano 
solar  gravity  dominated  the  nodal  motion.  For  long  periods  it  was  shown  that  the  orbit 
inclination  and  nodal  orientation  each  exhibited  an  oscillatory  behavior  having  the  same 
period.  The  findings  revealed  a  resonance  exists  between  the  inclination  and  the  nodal 
motion  [68]. 


5.3  The  Nodal  Motion 

To  study  the  long  term  nodal  motion  of  Sun  synchronous  orbits  it  is 
sufficient  for  this  investigation  to  determine  the  results  to  0(J2).  The  rate  of  change  jf 
the  node  from  (2.17)  is 


Q  =  "3n/2  cos  i  J2  (Re/p)2 


(5.1) 
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Let 

£2c  =  the  clock  angle,  the  angle  between  the  longitude  of  the  ascending  node 
and  the  mean  Sun  used  computing  the  local  time. 

0S  =  the  desired  precession  rate  of  the  longitude  of  the  ascending  node;  the 
same  rate  as  the  apparent  eastward  movement  of  the  Sun. 

Rg  =  mean  equatorial  radius  of  the  Earth 

Start  with  the  relationship 

fic  =  6-0s  (5.2) 

Differentiating  (5.2)  and  remembering  0S  =  constant 

Qc  =  fi  =  3n/2  sin  i  J2  (Re/p)2  di/dt  (5.3) 

Since  the  change  in  inclination  due  to  even  zonal  harmonics  is  zero,  the  next 
area  to  investigate  is  the  change  due  to  solar  gravitation.  The  mathematics  of  the 
gravitational  effects  of  the  Sun  and  Moon  on  near-Earth  satellites  was  first  developed 
by  Musen,  et  al.  in  1961  [70],  These  results  were  later  adapted  to  a  more  convenient 
form  which  allowed  the  disturbing  function  to  be  expressed  in  terms  of  osculating 
Keplerian  elements  for  use  in  the  equations  of  motion. 

The  complete  gravitational  disturbing  function  due  to  a  body  such  as  the  Sun 
or  Moon  orbiting  the  Earth  as  derived  by  Kaula  [71]  is 

R=^*?  7^1  5.p,h.qjkmoSF"mp(i)Fnrnh(i*)Hnpq(e)Gnhj(e*) 

x  cos[(n-2p)  co+  (n-2p+q)M  -  (n-2h)  to*  -  (n-2h+j)M* 

+  m(Q-ft*)]  (5.4) 
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where  F,  H,  and  G  are  polynomials.  For  tables  of  Fnmp/h  and  Gnhj  see  Kaula  [12J.  A 
table  for  Hnpq  is  contained  in  reference  72.  The  elements  of  the  disturbing  body  are 
designated  by  asterisks  (see  Figure  5.1).  Also, 

kQ  =  1  and  kj,,  =  2,  m  *  0 

Examining  the  trigonometric  expression,  the  terms  likely  to  be  significant  are 
those  of  long  period.  When  n  -  2p  +  q  =  0  the  periodics  in  M  are  eliminated. 


\^_^/Sun/Mooh 


Sun/Moon 


If  oc  =  the  argument  of  the  trigonometric  expression  in  (5.4)  with  the  above 
values  of  the  n,  p,  q,  and  j  integers  then 

«  =  [-(n  -  2h)cb*  -(n  -  2h)M*+  mQ]  (5.5) 

For  Sun  synchronism  when  the  Sun  is  the  disturbing  force 
(n  -  2h)  M*  =  mG 
Therefore  from  (5.5) 


n  -  2h  =  m 


For  long  term 


Since 


n  -  2h  +  j  *  0 
n  =  2  and  j  =  0,  h  *  1 
h  =  0  then  m  =  2 


Substituting  the  determined  values  of  the  integers  n,  p,  q,  j,  h,  and  m  into 


12a*3  F;221  (i)F22o(i*)H2io(e)G2oo(e*) 
x  cos  [2gj*-  2M*+  2(Q-n* )] 


From  the  tables  of  F,  G,  and  H 


F220  0*)  =  3(1+  COS  i*)2  1 4 
F22i  (i)  =  3/2  sin2i 
G200  (e*)  an^  ^210  (e)  ~  ^ 


Expressing  the  change  in  inclination  in  terms  of  the  disturbing  function 


[cos  i  ^ 


(pa) /2  (1- e2) /2  sin  i  dw  dQ 


A, 
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Substituting  (5.7)  into  (5.6)  and  taking  the  appropriate  partials 

&=0  (5.9) 

dco 


=  -3p*a2  sin2i  (i+  cos  i*)2sin  (-2(co*+M*+Q*-Q)) 
0“  16a*3 

(5.10) 


From  figure  5.1  Qc  =  to*  +  M*  +  Q*  -  Q 


(5.11) 


Substituting  (5.9),  (5.10),  and  (5.1 1)  into  (5.8)  the  change  in  inclination  due 
to  solar  attraction  is 


di  3  u*  2 

dt  =  "l6na*^  sin  1  (1+  cos  **)  sin(-2Qc) 


(5.12) 


From  IAU  (1976)  System  of  Astronomical  Constants  for  the  Sun 
i*  =  solar  obliquity  =  23.45° 
p*  =  1.32715  x  10nkm3/sec2 
a*  =  1.49467  x  108  km 

Equation  (5.12)  compares  exactly  with  the  analytically  averaged  results 
found  in  refeience  68  in  which  the  model  used  consisted  of  even  zonal  harmonics 
through  order  4  and  solar  gravitation. 

The  libration  period  of  the  node  can  now  be  found.  For  regions  where  drag 
is  negligible  and  assuming  that  sin2i  is  approximately  constant  substituting  (5.12)  into 
(5.3)  hc  =  Asin2Qc  (5.13) 

A  =  32a*}  J2  (^£)2(l+cosi*)2sin2i 
p  =  a(l-  e2) 


where 


Equation  (5.13)  is  in  the  form  of  the  familiar  pendulum  equation  and  can  be 
solved  analytically  using  elliptic  functions.  The  period  is  found  to  be 


T=  2k  [1+  !/4  sin2  £3C0  +  9/64  sin4  Qco  +  25/256  sin6  Qco  +  -  -]/(2c)1/2 

(5.14) 

where  Qqq  =  the  initial  clock  angle 

C=  "32a*5  3 2  (1+  cos  i*)2sin2  i 


The  actual  veracity  of  the  findings  in  this  section  can  be  demonstrated  by 
comparison  with  actual  1985  NORAD  data. 

Comparison  #1:  #1 1416  U-  NOAA6  Weather  Satellite 


*1  = 

006.91066734  day 

*2” 

364.51214982  day 

*1  = 

98.5704° 

h  = 

98.5375° 

el  = 

.0012005 

e2  = 

.0011415 

al  = 

7187.775  km 

*2  = 

7187.295  km 

^Cl= 

110.56° 

oC2= 

106.52° 

(1922.24  L  ascending  node) 

(1906.08  L) 

From  (5.12)  <^1=  -.085139  x  1  O'3  de8/day  di/dt=  -.070557  x  10"3  de8/day 

AiAvg  =  -.077848  x  10'3  **/day 

At  =  357.6014825  days 

<Ai,o„p  Avg  =  --02784'  vs.  AiacIua,  =  -.0329°) 

Comparison  #2:  #12553  -  NOAA7  Weather  Satellite 


<1  = 

020.461321 10  day 

h  ~ 

365.43820864  day 

h  = 

99.0262° 

x2  = 

99.0747° 

e,  = 

.001327 

e2  = 

.0012216 

ai  = 

7227.384  km 

a2  = 

7226.827  km 

QC1=  51.25° 

(1525  L  ascending  node) 


QC2=  60.25° 

(1600.54  L) 

From  (5.12)  &/<&=  .127187  x  10"3  ^S/day  <tydt=  .112291  x  10'3  ^S/day 

AiAvg=  .1 19739  x  10-3de8/day 
At  =  344.9768875  days 

<AicompAvg=-0413°  vs.  Aiactual=  .0485°) 

Table  5-1  CLOCK  ANGLE  LIBRATION  PERIOD 
(initial  orbital  elements  for  NOAA6  and  NOAA7  used) 


^co 

(deg) 

NOAA6  Period 
(years) 

NOAA7  Period 

(years) 

0 

24.070 

24.233 

15 

24.489 

24.655 

30 

25.823 

25.997 

45 

28.219 

28.409 

60 

31.479 

31.692 

75 

34.540 

34.774 

90 

35.823 

36.065 

To  study  the  long  term  effect  of  lunar  gravitation  on  the  change  in  inclination 
the  same  steps  are  taken  as  were  followed  for  the  examination  of  solar  gravitation. 
Using  identical  reasoning  n  =  2,  q  =  0,  p  =  1,  and  j  =  0.  The  argument  of  the 
trigonometric  expression  in  (5.4)  is 

<=*  =  [(n-2p)  (u+  (n-2p+q)  M-  (n-2h)  co*  -  (n-2h+j)  M*  +  m  (D-Q*)J 
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Using  these  values  for  n,  q,  p,  and  j  equation  (5.8)  becomes 


di  _E! _  f  t  (2-m)! 

a3 

x  H210(e)G2h0(e*) sin« 


(i)F2mh  C*) 
(5.16) 


where  «  =  [  (2h-2)(to*+M*)  +m  (Q  -  Q*)]  (5.17) 

Equation  (5. 16)  can  be  written  in  the  following  form 

S  sin  oe 
dt 

where  S  is  the  summation  over  the  possible  integer  values  of  m  and  h.  From  the  1984 
Astronomical  Almanac  the  following  is  data  for  the  Moon: 
e*  =  .054900489 
a*  =  3.84402  x  105  km 
p*  =  .0048994  x  106  km3/sec2 

i*  =  varies  over  an  18  year  period  from  18.32°  to  28.58°; 

average  value  will  be  used  in  this  study. 

Q*  =  period  of  18.6  years 
to*  =  period  of  8.9  years 
Q  =  period  of  1  year 
M*  =  period  of  27.21222  days 

Table  5-2  lists  the  values  of  S  for  the  respective  long  term  integers.  For  a 
complete  listing  of  values  of  S  containing  all  possible  integer  values  m,  p,  h,  and  j  see 
Appendix  C. 


Table  5-2  LUNAR  GRAVITATIONAL  EFFECT  ON 
NEAR-CIRCULAR  ORBITS 
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n  =  2,  p=  1,  j  =  0,  q  =  0 

m _ b _ 


o 

0 

1 

1 

2 

2 


0 

0 

0 

2 

0 

2 


S  fdeg/dav) 

0 

0 

.186543  x  10-4 
-.803530  x  10‘6 
.282902  x  10'3 
.524913  x  106 


The  only  value  that  has  an  effect  on  (5.16)  of  the  order  of  magnitude  of  that 
of  the  Sun's  gravitational  effect  is  S  =  .282902  x  10' 3  sin  «.  From  (5.17) 

sin  «=  =  sin  [-2  (co*+  M*+  Q*-  Q)]  (5.18) 

Since  M*  has  a  period  of  27.21222  days  this  is  the  controlling  orbital 
element  and  over  a  long  period  of  time  the  other  elements  have  little  effect  on  the  value 
of  the  trigonometric  expression.  Approximately  once  a  month  the  value  of  the 
trigonometric  term  will  have  cycled  through  one  period  and  the  effects  will  cancel. 
Therefore,  over  the  long  term  the  effects  of  lunar  gravitation  can  be  ignored  in 
comparison  to  the  effects  of  the  Sun. 
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5.4  Repeated  Ground  Track 

For  altimetry  applications  it  is  often  desirable  for  the  satellite  to  periodically 
repeat  its  ground  track.  This  was  accomplished  by  SEASAT  A  during  its  last  month  of 
operation  and  the  Earth  Observatory  Satellite  (EOS).  Additionally,  repeat  ground 
tracks  are  planned  for  many  altimeter  missions  now  under  consideration.  Ground  track 
repeatability  is  important  for  satellite  science  experiments  which  are  correlated  to 
ground  measurements  or  which  require  multiple  samplings  over  a  particular  point  on 
the  Earth.  Repeating  the  track  furnishes  a  predictable  pattern  of  coverage  and  permits 
direct  comparison  of  similar  data  taken  at  regular  intervals. 

Sun  synchronous  orbits  permit  the  repetition  to  occur  under  similar  lighting 
conditions.  The  interval  of  time  during  which  Earth  coverage  is  fully  completed  by  a 
satellite  is  referred  to  as  the  repeat  cycle ,  N,  which  is  measured  in  days.  The  number 
of  orbits  during  a  repeat  cycle  is  m,  an  integer.  The  value  of  N  is  constrained  to  be  an 
integer  if  the  initial  track  of  a  repeat  cycle  is  required  to  be  superimposed  on  the  initial 
track  of  a  previous  cycle  at  the  same  time  of  day.  Each  combination  of  N  and  m  is 
valid  if  it  satisfies  the  requirement  for  periodicity  and  similar  lighting  conditions, 
complies  with  altitude  restrictions  (imposed  by  launch  vehicle,  experiment 
instrumentation,  tracking,  orbit  adjustment,  etc.),  and  achieves  full  Earth  coverage. 
Full  Earth  coverage  will  result  if  the  sensor  swath  width  is  adequate  for  the  selected 
value  of  m. 

If  the  requirements  restrict  the  semimajor  axis  to  a  range  of  values  the  bounds 
on  the  mean  motion  and  period  of  the  orbit  are  determined.  For  this  study 

7115  km  <  a  <  7500  km.  Therefore, 

nmax=v/^"  =  5208.052°/day  =>  T  =  99.546  minutes 
nmin=  4812.223°/day  =*  T  =  107.734  minutes 
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This  in  turn  places  limits  on  the  number  of  revolutions  per  day. 

13.367  <  #  of  revs  <  14.467 

The  desired  range  of  semimajor  axis  translates  into  bounds  on  the  number  of 
revolutions  per  day. 

If  there  is  a  requirement  for  a  repeat  cycle  of  N  days  it  should  be  remembered 
that  an  integer  number  of  revolutions  are  necessary.  A  fraction  with  the  numerator  as 
the  number  of  revolutions  in  the  repeat  cycle  time  and  the  denominator  as  the  whole 
number  of  days  for  the  repeat  cycle  should  be  formed.  This  fraction  which  represents 
the  number  of  revolutions  per  day  must  be  within  the  determined  limits.  If  the  fraction 
is  reducible,  the  repeat  cycle  is  actually  less  than  the  desired  N  [73]. 


5.5  Nodal  Distance 

The  first  order  gravitational  term  (J2)  is  normally  the  most  significant 
perturbation  to  the  nodal  period  for  near-Earth  satellites.  For  orbital  motion  above  400 
km  (see  Tables  4-4  and  4-5)  the  effect  of  atmospheric  drag  is  much  smaller  than  that 
introduced  by  J2.  Therefore,  for  the  determination  of  the  distance  between  consecutive 
ascending  nodes  in  an  Earth-fixed  frame  only  the  effects  of  J2  will  be  considered  here. 

The  nodal  period,  PN,  is  the  time  necessary  for  the  argument  of  latitude  to  go 
through  2k  radians.  The  anomalistic  period,  PA,  is  the  time  to  go  from  perigee  to 
perigee.  Assuming  no  secular  change  in  a,  e,  and  i 


let 


PA  -  2k  (a%)>'2 


(5.19) 
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PN  can  be  obtained  from  the  proportion 


PN  =( 


2l L 


2rc  +  co  PA 


)Pa 


I 

(5.20)  ! 


where 


cb  =  change  in  argument  of  perigee 
PN  =  PA  in  the  absence  of  perturbations 


Also,  from  references  56  and  74 

DN  =  PN(C0e-a)  (5.21) 

where 

Dn  =  distance  between  consecutive  ascending  nodes; 
nodal  distance 

Q.  =  inertial  nodal  precession  rate 
coe  =  Earth  rotation  rate 


Substituting  (5.19)  and  (5.20)  in  (5.21)  and  simplifying 


Dn  = 


2n(jr) /2  (coe- a) 


(5.22) 


For  this  section,  let 

P  =  0-  e2)2 
f  =  (1-  5/4  sin2i  ) 

R  =  radius  of  the  Earth  at  the  equator 
C  =  J2  R2 


m,  N  are  defined  in  Section  5.4 


Using  equations  (2.17)  and  (2.18)  for  to  and  Q  to  0(J2)  and  making  the 
appropriate  substitutions,  the  nodal  distance  for  a  perturbed  orbit  is 


2  ( 7T  )  /2  p  coe  +  3  C  cos  i 

- 1(5'23) 


Rearranging  (5.23) 


2  Dn  a2p  -  4k  (a7/p)1/2  p  to  =  67t  C  cos  i  -  6DN  Cf  (5.24) 


Substituting  h2  =  a  in  (5.24) 


— P— h 7 . 2Dn  ph4  +  6C(DNf  -  rc  cos  i)  =  0  (5.25) 


One  way  to  achieve  the  desired  repeat  cycle  while  taking  into  account  the 
gravitational  effects  of  J2  is  to  determine  the  desired  nodal  distance  for  an  unperturbed 
satellite.  The  semimajor  axis  that  is  necessary  to  achieve  the  same  nodal  distance  under 
the  influence  of  the  gravitational  perturbation  can  then  be  computed. 

The  steps  to  be  taken  are: 

1.  Compute  the  nodal  distance  for  unperturbed  motion  as  follows 
PA  (mins)  =  (24  x  60)N/m 
The  semimajor  axis  for  unperturbed  motion  is 


AD-A17*  679  FROZEN  ORBITS-NEAR  CONSTANT  OR  BENEFICIALLV  VARYING 
ORBITAL  PARAAETERS(U)  AIR  FORCE  INST  OF  TECH 
HRIGHT-PATTERSON  AFB  OH  R  C  HURROM  19  AAV  86 
UNCLASSIFIED  AFIT/CI/NR-86-89D  F/B  2/2 


2.  Given  an  eccentricity  and  inclination  angle  and  using  the  required  DN  in 
(5.25)  solve  the  seventh  order  polynomial  to  determine  the  required  semimajor  axis. 

3.  If  a  Sun  synchronous  orbit  is  desired  an  iterative  process  will  be 
necessary.  For  step  #2  an  approximate  inclination  angle  can  initially  be  used  in  (5.25) 
to  determine  a  semimajor  axis.  The  inclination  angle  is  then  computed  using 

i  =  cos'1  [-4.7733  x  10'15  (1-  e2)2  a7/2]  (5.26) 

(assumes  =  .9856°/day) 

This  inclination  angle  is  then  used  in  (5.25)  and  the  process  is  repeated  until  the 
semimajor  axis  and  inclination  deviations  are  within  desired  tolerances. 

4.  The  computed  semimajor  axis  is  used  in  (5.19)  to  determine  the  new 
anomalistic  period. 


5.6  16-Dav  Repeat  Cycle 

For  a  16-day  repeat  cycle  the  fraction  discussed  in  section  5.4  would  have  a 
denominator  of  16.  To  satisfy  the  orbital  semimajor  axis  criteria  the  number  of 
revolutions  per  day  must  be  between  13.367  and  14.467. 

21Vl6  <  13.367  revs/day  <  2 1 4/ 1 6 
23 1  / 1 6  <  14.467  revs/day  <  232/ 16 

Since  213  revolutions  in  the  16  day  period  is  less  than  13.367  revs/day  it  is 
not  acceptable.  Additionally,  232  revolutions  in  16  days  would  necessitate  a  lower 
semimajor  axis  than  is  allowable.  Therefore, 


214  <  m  <  231 


Table  5-3  shows  each  fraction  with  a  numerator  between  214  and  231  and  a 
denominator  of  16.  If  the  fraction  is  reducible  the  orbit  does  not  qualify  as  a  16  day 
repeat  cycle.  Steps  1,  2,  and  4  in  the  previous  section  are  used  for  Table  5-3.  The 
semimajor  axes  listed  under  the  perturbed  column  are  those  required  to  yield  the 
desired  nodal  distance  for  e  =  .002  and  i  =  55.0°. 

Steps  1  through  4  are  used  to  generate  Table  5-4.  The  results  show  that  for  a 
Sun  synchronous  orbit  influenced  by  gravitational  perturbation  the  semimajor  axis 
required  to  produce  the  desired  nodal  distance  is  approximately  10  km  higher  than  the 
unperturbed  semimajor  axis. 
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Table  5-3  CHARACTERISTICS  OF  A  16-DAY  REPEAT  CYCLE 
(  )  =  valid  16  day  cycle  <  >  =  semimajor  axis 


e=.002  i=55.0° 

#  Revs  Repeat  cycle 

/1 6  days  Fraction  (days) 


UNPERTURBED 
i  i 

Orbital  Period 

Ht.(km)  (mins) 


PERTURBED 
i  A 

Orbital  Period 

Ht.(km)  (mins) 


13  7/16 

(16) 

1095.331 

<7473.494> 

107.163 

2990.509 

1051.501 
<7429. 664> 

106.222 

13  8/16=27/2 

2 

139/16 

(16) 

1049.325 

<7427.488> 

106.175 

2962.938 

1004.779 

<7382.942> 

105.221 

13  10/16=109/8 

8 

13  1 1/16 

(16) 

1004.018 
<7382.1 8 1  > 

105.205 

2935.869 

958.750 

<7336.913> 

104.239 

13  12/16=55/4 

4 

13  13/16 

(16) 

959.417 
<7337. 580> 

104.253 

2909.302 

913.422 
<7291. 585> 

103.274 

13  14/16=nl/8 

8 

13  15/16 

(16) 

915.479 
<7293. 642> 

103.318 

2883.210 

868.753 
<7246.9 16> 

102.327 

14  =  14/1 

1 

14  1/16 

(16) 

872.210 
<7250. 373> 

102.400 

2857.592 

824.749 
<7202. 912> 

101.396 

14  2/16  1 13/8 

8 

14  3. 16 

(16) 

829.570 
<7207. 733> 

101.498 

2832.410 

781.369 
<7159. 532> 

100.482 

14  4/16=-*7/4 

4 

14  5/16 

(16) 

787.500 
<7 165. 663  > 

100.611 

2807.668 

738.567 
<71 16.730> 

99.582 

14  6  16= 1 15  8 

8 

14  7  16 

(16) 

746.100 

<7124.263> 

99.740 

2783.361 

696.398 
<7074.561  > 

98.69S 

Table  5-4  CHARACTERISTICS  OF  A  SUN  SYNCHRONOUS 
16-DAY  REPEAT  ORBIT 
<  >  =  semimajor  axis 
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e=.002  UNPERTURBED  PERTURBED 


ft  Revs  Orbital 

16  days  Ht.(kirO 

Period 

(mins) 

Inclination 

(degs) 

dn 

(km) 

Orbital 

Ht.(km) 

Period 

(mins) 

Inclination 

(degs) 

215 

1095.331 

<7473.494> 

107.163 

99.918 

2990.509 

1106.610 
<7484. 773> 

107.406 

99.971 

217 

1049.325 

<7427.488> 

106.175 

99.704 

2962.938 

1060.486 
<7438. 649> 

106.414 

99.756 

219 

1004.018 
<7382. 1 8 1  > 

105.205 

99.497 

2935.869 

1015.062 

<7393.225> 

105.441 

99.547 

T 

959.417 
<7337. 580> 

104.253 

99.296 

2909.302 

970.344 
<7348. 507> 

104.486 

99.345 

223 

915.479 
<7293. 642> 

103.318 

99.100 

2883.210 

926.294 
<7304.4 57 > 

103.548 

99.148 

225 

872.210 

<7250.373> 

102.400 

98.912 

2857.592 

882.915 
<7261. 078> 

102.627 

98.958 

227 

829.570 
<7207. 733> 

101.498 

98.728 

2832.421 

840.166 
<7218. 329> 

101.722 

98.773 

229 

787.500 

<7165.663> 

100.611 

98.550 

2807.668 

798.004 
<7 1 76. 1 67> 

100.832 

98.594 

231 

746.100 
<71 24.263  > 

99.740 

98.377 

2783.361 

756.481 

<7134.644> 

99.958 

98.420 

CHAPTER  VI 


GEOSYNCHRONOUS  ORBITS 


6.1  Geosynchronous  Orbits 

Geosynchronous  satellites  orbit  approximately  36,000  km  above  the  surface 
of  the  Earth  at  a  speed  of  3.07  km/sec.  The  term  geosynchronous  refers  to  an  orbit  that 
has  a  period  equal  to  the  period  of  rotation  of  the  Earth  relative  to  an  inertial  system. 
Consequently,  the  satellite  is  at  rest  with  respect  to  the  rotating  Earth.  This  type  of 
orbit  is  mainly  used  for  communications  missions,  but  has  also  been  used  by  some 
Earth-observation  missions  and  scientific  missions. 

A  perfectly  geosynchronous  orbit  would  only  be  achiev  able  if  the  Earth  were 
perfectly  symmetrical  and  no  other  forces  were  acting  on  the  satellite  except  the  central 
gravity  attraction  of  the  Earth.  However,  additional  forces  acting  on  the  satellite  do 
change  the  shape  of  the  orbit,  the  orientation  of  the  orbital  plane,  and  the  longitude. 
These  changes  can  only  be  compensated  for  by  active  orbit  control. 

As  of  January  1,  1985,  there  were  approximately  185  satellites  that  are 
considered  near-stationary.  These  satellites  are  near-circular  (e  <  0.1),  near-equatorial 
(i  <  1(F),  and  have  an  orbital  period  close  to  one  mean  sidereal  day  (0.9  rev/day  <  n  < 

1 . 1  rev/day)  |3). 

6.2  Fundamentals 

With  the  origin  of  the  coordinate  system  at  the  center  of  mass  of  the  Earth  the 


m 


gravitational  potential  (2.1 1)  can  be  expanded  to  include  tesseral  harmonics  which  are 
due  to  longitudinal  variations  in  the  shape  of  the  Earth. 


u  -  4  !  l-f  Jn/f-' )"  Wn  4>)+|X  J Pnm(sin  40 

x  cos  m(  0  -6nm)]  (6.1) 


The  double  summation  includes  the  longitude-dependent  terms:  r  is  the 
geocentric  distance,  0  is  the  geographic  longitude,  4>  is  the  latitude,  Re  is  the  mean 
equatorial  radius  of  the  Earth,  0nm  is  the  longitude  in  the  direction  of  the  principal  axis 
of  symmetry  of  the  Earth's  distribution  accounted  for  by  the  nm  harmonic,  Pnm  is  the 
associated  Legendre  polynomial,  and  Jnm  is  the  numerical  coefficient  characterizing  the 
Earth's  mass  distribution  [4J. 

Using  a  spherical  coordinate  system  (figure  6.1)  that  is  rotating  at  the  same 
rate  as  the  angular  velocity  of  the  Earth  (to  )  the  force  field  F  of  the  Earth  on  a  point 


mass  is: 


c  c  A  c  A  c  A 

F  =  Fr  er+ F.e.  +  F„  e( 


/\  /\  /\ 

where  er ,  and  e0  are  unit  vectors  in  their  respective  directions. 


p,= 


(6.3a) 


(6.3b) 


r  _  m  d  U 
e  r  cos  <>  cFe" 


(6.3c) 


I 


EARTH 


nm  axis  of  mass  symmetry 
v,  of  the  Earth 


Vernal  Equinox 


Figure  6. 1  THE  COORDIN ATE  SYSTEM 


^■nm  \, 


.  Longitude  at  time,  t 


\  Vf-  nm  axis  of  Earth  symmetry 

/  at  time,  t 

Vernal  Equinox  y  T  <—  Greenwich  longitude  at  time,  t 

X  Greenwich  longitude  at  time  zero 


Figure  6.2  THE  EARTH  ELLIPSOID'S  EQUATOR 
WITH  LONGITUDE  REFERENCES 
(North  Pole  looking  south) 
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Figure  6.2  shows  the  ellipticity  of  the  Earth's  equator.  X  is  the  geographic 
longitude  of  the  satellite  while  Xnm  is  the  geographic  longitude  of  the  principal  nm  axis 


of  Earth  symmetry. 

It's  seen  that 

0nm  =  eG0  +  “el  +  Km 

(6.4) 

and 

x  -  jl  =  0  -  enm 

nm  nm 

(6.5) 

6-3  Equatorial  Near-Circular  Orbits 

The  short-period  variations  in  r,  0,  and  0  are  small  for  a  geosynchronous 
orbit.  On  the  equatorial  plane  there  is  a  contribution  to  F0  from  the  Jn  terms  with  n  odd 
and  the  Jnm  terms  with  (n-m)  odd.  The  effect  is  to  displace  the  orbital  plane  slightly 
from  the  equatorial  plane.  The  displacement  due  to  J3  is  approximately  50  cm 
southwards  from  the  equator.  From  this  fact,  F^  can  be  ignored  without  introducing 
any  meaningful  error  [41  ]. 

In  addition,  FR  can  be  ignored  for  this  study.  The  effect  of  the  zonal 
harmonics  can  be  observed  using  (6.3a)  The  semimajor  axis  corresponding  to 
synchronous  motion  is  reduced  by  only  .5  km  if  J2  is  ignored  [41  ]. 

The  tesseral  harmonics  primarily  affect  the  tangential  force  F0.  Although 
these  forces  are  small,  their  long-term  effects  can  be  large  when  the  motion  of  the 
satellite  is  commensurable  with  the  Earth's  rotation  rate. 

Employing  an  approach  used  by  Allan  (1963)  the  longitude  acceleration  can 
be  determined.  The  total  unperturbed  energy  of  the  satellite  in  a  spherical  Earth  gravity, 
field  is 


therefore 


dE  =  M/2a2  da 


(6.6) 
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The  energy  change  caused  by  the  force  acting  on  the  satellite  each  day  is 
approximately  constant 


f 2- 

dE  =  J  F  •  ade  =  2na  F  (6.7) 


For  an  orbit  that  is  near-circular  the  averaged  force  (F)  per  orbit  is  mainly  the 
tangential  perturbation  force  Fe. 

Substituting  Fe  =  F  in  (6.7)  and  then  equating  the  result  with  (6.6) 
u/2a2  da  =  2na  Fe  (6.8) 

Solving  (6.8)  for  da,  recalling  p  =  n2a3,  and  using  an  orbital  period  of  1 
sidereal  day 

da  =  2/n  Fe  (6.9) 

Starting  with  p  =  n2a?  and  using  implicit  differentiation 

a  =  '2"a/3n  (6.10) 

Differentiating  9  =  n  -  nu  (6. 1 1 ) 

where  n  is  the  true  mean  motion  and  no  is  the  mean  angular  velocity  of  the  Earth 

0  =  n  (6. 12) 


Substituting  (6. 12)  in  (6. 10)  and  combining  with  (6.9) 
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0  = '3/a  F0  (6.13) 

It  is  seen  from  this  equation  that  the  apparent  long-period  acceleration  0  is 
always  opposite  to  the  longitudinal  force.  Equation  (6.9)  shows  that  if  the  force  is  in 
the  direction  of  motion  (i.e.  +)  then  the  semimajor  axis  is  increasing  as  is  the  energy. 
At  the  same  time  the  mean  motion  and  0  are  decreasing. 

Since  F0  gives  the  long-period  change  in  the  semimajor  axis,  it  also  changes 
the  period  of  the  satellite.  From  Kepler's  third  law,  the  period  of  a  24-ho’ir  orbit  is 

T  =  2Ji(a3/|i)1/2  (6.14) 

Therefore 

dT  =  3>n  (a/|i)1/2  da  (6.15) 


Substituting  (6.9)  for  da  in  (6.15)  the  change  in  period  each  day  of  a 
24-hour  near-circular  orbit  is  given  by 

dT=  12*-’ (aW:  F„  (6.16) 

If  the  present  investigation  is  restricted  to  the  dominant  J22  term  m  =  n  =  2 
in  (6.1 )  and  (6.3c)  then  substituting  into  (6.13)  with  <J>  =  0,  r  =  a,  and  n  =  2n 

0  =  A22  sin  2  (0  -  022)  racVsid.day2  (6.17) 
where  A22  = -72’t2  J22  (^c/a)2 

The  coefficient  A^  is  always  positive  if  all  the  tesseral  Jnn  terms  in  the 
gravity  potential  are  arbitrarily  assigned  as  negative  numbers.  Equation  (6.17) 


represents  the  drift  acceleration  as  a  function  of  the  longitude  in  a  second-order  gravity 
field  for  an  equatorial  synchronous  satellite  where  F0  is  constant  at  every  point  in  the 
orbit.  Satellite  data  has  shown  that  for  eccentricities  as  high  as  .0012  the  equation  is 
extremely  accurate.  If  9  is  taken  as  the  mean  longitude  of  the  ascending  and 
descending  nodes  (6.17),  with  reasonable  accuracy,  represents  the  24-hour  drift 
regime  for  eccentricities  as  high  as  0.3  [49]. 

The  results  of  this  approach  compare  exactly  with  those  of  the  method 
completed  by  Kaula  [  12)  in  which  the  following  resonating  disturbing  function  was 
used 

R=(n^m)  QnmCOSm(0)+  M  *  Q'  '  9  C0 )  (618) 

'even 

where  Qrm  is  a  polymonial  function  of  inclination,  eccentricity,  semimajor 
axis,  and  respective  spherical  harmonics  (see  page  50,  reference  12).  R  can  be 
exercised  in  the  Lagrange  planetary  equation 

da/dt  2  na  dR/5M 

The  change  in  semimajor  axis  will  cause  an  acceleration  in  the  satellite's 
longitude  which  is  defined  by  K  =  to  -  \1  +  Q  ->  0G  (6.19) 

Using  Kepler's  law 

M  =  ndt  =>  M  =  dn/dt 

Therefore 


/,  - 


M  =■  cjn/0a  da/dt  =  2/na  dn/da  fjR/r)\1  =  -?/a:  3r/0M 


(6.20) 


Taking  the  partial  of  R  (6.18)  with  respect  to  M,  using  (6.19)  and 
substituting  in  (6.20)  the  result  is 

mQira  sinm(X-Xnn,) 

M  v"  “''even 

when  n  =  m  =  2 

X  =  Q22  sin  2  (  X- Xnm)  (6.21) 

Q22  =  3n2  a2  J22  (Re/a)2  @  i  =  0° 

letting  n2  =  4rc2 

Q22  =  - 1 2xc2  a2  J22  (Rc/a)2  (6.22) 

Substituting  (6.21 )  and  (6.5)  in  (6.20) 

e  =  X  =  -72k2  j22  (Rc/a)2  sin  2  (0  -  0nm)  (6.23) 

which  is  equivalent  to  (6.17). 

I  "sing  (6. 10)  and  (6. 1 1 ) 

0  =  "3n/2a  da  (6.24) 


Solving  (6.8)  for  da  and  substituting  in  (6.24)  with  n  -  2 n 


I 
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A  _  -127r2a2  p 
0-  4  he 


Substituting  for  F0  to  second-order  with  <J)  =  0° 

0  =  -24tt2  J22  (^e/a)2  s*n  ^  (0  • 

0  =  ^22/3  sin  2  (0  -  022)  racVsid.  day  (6.25) 

This  represents  the  apparent  net  longitudinal  drift  rate  of  the  24-hour 
equatorial  near-circular  orbit's  ground  track  with  respect  to  the  surface  of  the  Earth. 

6.4  Stable  and  Unstable  Equilibrium  Points 

The  actual  tangential  force  on  the  satellite  from  (6.13)  and  (6.17)  is 

Fe  =  'a/3  0  =  'a/3  A22  sin  2  (0-0,,)  (6.26) 

It  is  seen  from  this  equation  that  the  force  vanishes  when  0  -  022  =  nl 2  1  where  1  =  0, 
1,  2,  and  3.  These  values  correspond  to  the  major  and  minor  axes  of  the  equatorial 
ellipse.  The  force  is  directed  towards  the  nearest  end  of  the  major  axis  with  the 
maximum  value  at  the  45°  points  between  the  axes. 


i  .  *  <  ‘j»  *  *  •  '  >  *  <  fc  •  "  ■»' 


. .  .  ,  _  .  .  .  ,  .  %  ,  ....  „  «,  „  .  ,  ,  .  < 


Figure  6.3  THE  EQUILIBRIUM  POINTS 

Blitzer,  et  a!.,  provided  a  mathematical  derivation  demonstrating  that  the  two 
points  on  the  minor  axis  are  positions  of  stable  equilibrium  while  those  on  the  major 
axis  are  unstable  points  [39],  Later,  Blitzer  showed  that  when  the  entire  spectrum  of 
tesseral  harmonics  is  included  the  equilibrium  points  were  no  longer  symmetrically 
located  on  the  extensions  of  the  principal  axes  of  the  equatorial  ellipse.  The  radial  and 
latitudinal  displacements  from  the  symmetrical  positions  are  negligible  while  the 


longitudinal  deviations  can  vary  by  3°  or  more.  The  dominance  of  the  J22  term  still 
limits  the  number  of  equilibrium  positions  to  four  [44],  The  stable  longitudes  are  near 
75'E  and  105°W  whereas  the  unstable  longitudes  are  approximately  165°E  and  15°W 
[41,75], 

Since  the  stable  equilibrium  points  exist  along  the  minor  axes,  it  is 
convenient  to  change  the  origin  (see  figure  6.3).  Thus, 

V  =  0  -  022  ±  nl2 

Equation  (6.17)  becomes 

\|/  =  -A22  sin  2  \\i  (6.27) 

Multiplying  (6.27)  by  x\f  leads  to  a  first  integral  which  can  be  written  as  [41  ] 
V2  -  2A22  cos2  y  =  C  (6.28) 

C  is  a  constant  of  the  motion  which  is  determined  by  the  initial  conditions. 
From  (6.28)  it  follows  that  C  >  -2A-,-,  cos2  yo.  If  the  initial  conditions  are  such  that  C 
is  small,  then  \\i  cannot  go  through  a  full  cycle.  If  C  is  such  that  there  is  a  maximum 
departure  of  ym  which  is  captured  between  'K/2  and  n/2  then  \|/  =  0. 

Therefore,  from  (6.28) 

C  =  -2A22  cos2  \j/m 

and 

cos  vj/m  =  (C/-2A22)1/2  (6.29) 

If  the  maximum  value  achieves  is  n/2  then  C  must  equal  zero.  If  vj r  is 
less  than  n/2  then  C  <  0. 

Referring  back  to  (6.28)  with  C  <  0 


Vo2  -  2A22  cos  2  V0  -  0 

Therefore, 

I V0 1  <  I  (2A22)1/2  cos  v0  | 


100 


(6.30) 


If  this  condition  is  satisfied  the  satellite  will  oscillate  in  longitude  about  the 
stable  equilibrium  point.  If  not  satisfied,  the  satellite  circulates  continuously  around  the 
Earth's  longitudes.  In  addition,  the  maximum  angle,  ±Vm,  to  which  the  satellite  is 
confined  in  its  oscillation  is  found  from 

V02  '  2A22  cos  2  v0  =  -2A22  cos2  vm 

Therefore, 


Using  (6.28)  with  Vm  =  0 


V2  =  2A22  (cos2v  -  cos2Vm)  (6.32) 


To  find  the  period  of  the  oscillation  proceed  in  the  following  manner: 
substitute  cos 2\\i  =  1-  sin2v  in  (6.32)  and  solve  for  v 

V  =  (2A22),/2  (sin2vm  -  sin2v)1/2 

let  k2  =  Vsin2vm 

therefore, 

V  =  ^V/dt  =  (2A22),/2  ('/j^  -  sin2v)1/2 


(6.33) 


m 


Table  6-1  PERIOD  OF  OSCILLATION  ABOUT 
THE  STABLE  EQUILIBRIUM  POINT 


J22  =  -1.7  x  10" 
a  =  42,163  km 
i  =  0° 


Vo=0  =>  V0  =  ¥m  (see  6.31) 
y  =  sin1  (k1/2  sin  \j/m) 

A22  =  2.764  x  10'5  rad/sid. day 


Y  (deg 


4.627 

1291.041 

5.791 

919.916 

5.000 

760.833 

3.297 

671.023 

1.073 

614.676 

8.529 

578.108 

5.784 

554.991 

2.920 

542.117 

1013.21 

1353.40 


yfFiw 


This  type  of  orbit  has  the  following  particular  characteristics: 

(1)  The  satellite  crosses  the  equator  twice  a  day 
-  at  the  ascending  and  descending  nodes. 

(2)  Excluding  perturbations,  the  ground  track  is  repeated. 

(3)  The  ground  track  is  a  closed  curve. 

For  an  inclined  24-hour  satellite  the  energy  changing  force  on  the  satellite  is 

no  longer  constant  over  a  single  orbit.  A  satellite  in  a  near-circular  orbit  now  describes 
a  closed  figure  eight  path  over  the  Earth's  surface  centered  on  the  equator. 
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Figure  6.4  SATELLITE  GROUND  TRACK 


At  each  point  on  the  path  the  tangential  force  varies,  as  both  the  longitude  and 
latitude  change  during  the  daily  excursion.  The  force  along  the  track  is  now  composed 
of  contributions  from  both  latitude  and  longitude  gravity  perturbation  forces.  The  zonal 
gravity  forces  have  no  net  daily  energy  effect.  For  the  inclined  orbit  the  drift  equations 
of  motion  will  contain  many  small,  short-period  alterations  introducing  difficult 
non-linearities.  To  arrive  at  a  theory  for  long-term  drift  motion  Wagner  [49]  used  a 
first  approximation  perturbation  technique  which  smooths  out  the  short-period  effects 
and  averages  over  the  24-hour  period. 

Using  the  tangential  perturbing  force  (FT)  in  place  of  Fe  in  (6.13) 

0  =  '3/a  Ft  (6.38) 

where  FT  =  F0  cos  +  Ffi  sin  « 

=  angle  between  a  meridian  plane  and  the  orbital  plane 


and  proceeding  on  the  basis  of  the  orbit-averaging  theory  established  by  Wagner  in 
reference  49  an  inclination  factor  can  be  determined.  This  factor  is  applied  to  the  zero 
inclination  regime  to  get  the  proper  acceleration  for  the  inclined  satellite.  The  simple 
factor  which  is  independent  of  the  longitude  is  used  to  modify  the  equatorial  regime. 

The  orbit-averaging  modification  gives  the  complete  drift  through  m  =  n  =  4 


4  n 


X  =  -1271  V  V  C  F(i) 

^  nm  nm 


n=2  m=l 


j 


nm 


(*i.c  \  sin  m  (X  -  X  ) 

t  r  >  nm 


rad 

/sid.  day  ^ 


as 


C22  -  6 

NJ 

II 

( 1  +  cos  i ) 1 

4 

„  -3 

f(i)31  = 

2 

(1+cosi)  5  sin  i  (1+3  cos  1) 

C31  ~  2 

2  8 

C33  -  45 

F(i)33  = 

(1+  cos  i)3 

8 

C42  =  '15 

F(i)42  = 

2  2 

(1+ cos  i)  7  sin  i  cos  i  (1+  cos  1) 

4  4 

C44  =  420 

F(j>44  - 

(1+  cos  i)4 

16 

The  contribution  to  the  longitude  acceleration  of  the  inclined  24-hour 
near-circular  orbit  due  to  the  32, 41,  and  43  harmonics  is  zero. 

As  was  demonstrated  in  section  6.4  (see  Table  6-1)  the  amplitude  of  the 
oscillation  depends  on  the  initial  longitude.  Figure  6.5  shows  the  amplitudes  of  the 
oscillations  for  orbits  with  an  initial  longitude  of  75°  from  the  stable  point  on  the  minor 
axis  (75°E).  The  effect  of  inclination  is  to  reduce  the  rate  of  change  of  the  mean 
longitude  and  as  a  result,  increase  the  period  of  the  oscillation. 

1QA  —  Initial  Lomtitude  =  0  dci;. 


Longitude 


1  2  Time  (Years)  3 

Figure  6.5  MEAN  LONGITUDE  VS.  TIME  [75] 


Because  of  the  increasing  number  of  applications  of  the  geosynchronous 
orbit,  available  space  in  the  orbital  arc  is  being  depleted.  The  NASA  Space 
Transportation  System  is  expected  to  launch  approximately  160  geosynchronous 
satellites  between  1980  and  1990  [76].  These  satellites  need  to  have  the  proper 
separation  to  avoid  physical  and/or  radio  frequency  interference.  The  geosynchronous 
arc  of  interest  to  the  domestic  U.S.  is  55° W  longitude  and  135°W  longitude. 
Extensive  future  traffic  and  the  large  amounts  of  arc  needed  for  certain  projects,  such 
as  space  manufacturing,  will  promote  the  use  of  the  inclined  orbit.  In  addition,  inclined 
geosynchronous  orbits  are  useful  for  applications  where  there  must  be  extended 
viewing  time  [75]. 

An  enormous  advantage  of  the  inclined  geosynchronous  orbit  is  it  allows 
several  satellites  to  have  the  same  figure  eight  ground  track.  They  will  cross  the 
equator  at  the  same  longitude,  called  the  gateway.  The  gateway  may  require  only  a  few 
degrees  of  the  geosynchronous  arc,  while  accommodating  up  to  10-12  satellites  [75]. 

6.6  Velocity  Requirements  for  Orbit  Maintenance 

To  maintain  a  satellite  in  a  synchronous  position  other  than  the  two  stable 
positions  it  is  necessary  to  do  so  \sith  propulsion.  To  determine  the  AV  requned  to 
offset  the  orbit  perturbation  consideration  will  be  given  only  to  in-plane  corrections. 
The  following  assumptions  are  made  for  this  study:  (1)  the  orbit  is  near-circular  so 
that  AV  in  the  radial  direction  is  negligible  as  compared  to  the  tangential  component, 
and  (2)  there  is  no  coupling  between  the  in-plane  and  out-of-plane  perturbations. 

The  orbit-averaged  longitude  drift  of  an  inclined  near-circular  orbit  to  0(Jt,) 
is  obtained  from  (6.39) 


X  =  D-,.,  sin  2  (X  -  X-,-,)  racVsid.day-  (6.40) 


where 
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D22  —  -72n2  (^e/r)2  J22  F(i)22 

Following  the  approach  used  in  references  49  and  75  the  above  equation  can 
be  linearized  in  terms  of  the  variation  A X  defined  by 

X  =  XS  +  AX  (6.41) 

where  Xs  is  the  desired  mean  longitude  of  the  ascending  node.  Substituting  (6.41)  into 
(6.40)  and  assuming  AX  is  small  the  result  is 


AX  -  [2D22  cos  2(XS  -  Xj-,)]  AX  =  D22  sin  2(XS  -  X2,) 


(6.42) 


The  complete  solution  to  (6.42)  is 

AX  =  A  cos  wt  +  B  sin  wt  -  '/2  tan  2(XS  -  X22)  (6.43) 

where  A  and  B  are  constants  and  w  =  [-2D22  cos  2(XS  -  X22)]'2. 

To  find  the  constants  of  integration  it  is  assumed  that  there  is  no  orbit 
injection  error  (i.e.  AX  =  0  and  AX  =  0  at  t  =  0).  The  solution  (6.43)  becomes 

AX  =  V2  tan  2(XS  -  X12)[  cos  wt  -  1]  (6.44) 

Equation  (6.44)  can  be  simplified  by  expanding  cos  wt  and  inserting  the 
expression  for  w.  The  result  is 


•»*  1  ^ „  •"»  ^ *  *  a  *  *•*  «  1 *  »  ^ ~  a  "■  ft  •  •  ”  »  * 

^  ^  1*1  il  *L*  r  *  - 


AA  = 


(6.45) 


D,2  t  sin  2(As  -  A22) 

2 

In  reference  49  it  is  shown  that  the  rate  of  change  of  the  near-circular  radial 
distance  satisfies  the  equation 

r  =  B22  sin  2(A  -  A22)  ‘“8*  units/sid.day  (6.46) 

where 

B22  =  24ttr  (^c/j.)2  J22  F(i)22 

Equation  (6.46)  can  be  linearized  in  terms  of  the  variation  Ar  defined  by 
r  =  rs  +  Ar  (6.47) 

Substituting  into  (6.46)  the  expressions  for  Ar  and  A A  and  assuming  AA  to  be  small 

Ar  =  B22  sin  2(AS  -  A^)  +  2B22  cos  2(AS  -  A22)  AA  (6.48) 

Substituting  (6.45)  into  (6.48)  and  solving  the  differential  equation 

Ar  =  B22  sin  2(AS  -  A22)  [  1+  */3  D22  t2  cos  2(AS  -  A22)]  (6.49) 

(the  constant  of  integration  is  zero  since  Ar  =  0  at  t  =  0). 

To  determine  the  AV  required  for  station  keeping  (6.45)  is  used  to  compute 
the  time  of  drift.  At  that  time  it  will  have  a  radial  distance  of  rs  +  Ar.  For  a 
near-circular  orbit  the  velocity  relative  to  inertial  space  is  given  by 


V.  =  (rc  +  Ar)  cm  -  r  v 


(6.50) 


where 


coe  =  Earth  rotation  rate 
v  =  rate  of  change  in  true  anomaly 

To  return  the  satellite  to  the  original  value  of  Xs  it  will  be  necessary  to 
transfer  a  near-circular  orbit  radius  rs  -  Ar  with  a  velocity  V2  relative  to  the  inertial 
space  given  by 

V2  =  (rs  -  Ar)  coe  +  rs  v  (6.5 1 ) 

Minimum  energy  expenditure  would  be  achieved  by  a  Hohmann  transfer 
with  rs  +  Ar  being  the  apogee  distance  and  rs  -  Ar  the  perigee  distance.  The  transfer 
can  be  accomplished  by  applying  an  impulse  AVj  at  time  t  . 

r  rs '  Ar  1/2  i 

AVi =  vi  1<~ — >  •'] 

S 

Since  Ar/rs  «  1  =*  AV,  =  ‘ '/2  V,  Ar/rs  (6.52) 

At  time  t0  plus  one  half  of  a  sidereal  day  a  second  impulse  AV-,  is  applied  to  slow  the 
vehicle  dow'n  to  the  desired  velocity  V,. 

r  +  Ar  ]n 

AV2  =  V2  [  1-  (-5- - )  -  ]  =  J/2  V2  Ar/r^  (6.53) 

S 

The  vehicle  will  move  back  towards  the  initial  position  X  as  Ar  increases  to  zero. 

Substituting  (6.50)  into  (6.52)  and  (6.5 1 )  into  (6.53)  and  neglecting  second 
and  higher  order  terms  of  the  small  variations  the  following  results 


Therefore,  the  required  AV  for  station  keeping  is 

AV  =  |  AV,  |  +  |  AV2  |  =  |  Arcoe  |  (6.54) 


The  period  for  the  complete  cycle  can  be  determined  from  (6.45) 


T  =  2  [ 


2AX 


D22  sin  2(X$-  X22) 


1/2 


(6.55) 


Using  (6.49),  while  neglecting  the  second  term  (due  to  order  of  magnitude), 
(6.45),  (6.54)  and  (6.55)  the  velocity  change  per  unit  of  operating  time  is  therefore 


to,  sin  2(>.s  -  X22) 


AV  we  a21 

~T 


(6.56) 


Figure  6.6  shows  the  in-plane  velocity  corrections  necessary  to  compensate 


for  longitude  drift. 


Longitude  Relative  to  Stable  Point  (Degrees)  (\^  -  X.,,,) 

Ficure  6  6  IN-PLANE  VELOCITY  CORREC  1  ION 
FOR  LONGITUDE  DRIFT 


This  section  is  not  intended  to  be  a  comprehensive  study  of  the  effects  of 
solar  radiation  but  rather  an  introduction  to  the  complexity  of  the  topic.  Although  the 
non-gravitational  perturbation  is  usually  small,  for  geosynchronous  satellites  with 
large  solar  collectors  it  could  be  substantial.  For  the  Solar  Power  Satellite  (SPS),  the 
perturbation  due  to  solar  radiation  was  of  the  order  of  magnitude  of  gravitational  terms 
due  to  the  Sun,  the  Moon,  and  J2.  The  major  effect  of  this  perturbation  is  oscillation  of 
the  eccentricity  [37,  69,  77], 


The  perturbation  equation  for  eccentricity  in  the  Gaussian  form  is 


de  (l-e*- 


ui  na 


2  cos  v  +  e(  1  +  cos  v) 


1  +  e  cos  v 


where  fj  and  f2  are  the  radial  and  transverse  components  of  the  disturbing 


acceleration. 


For  near-circular  orbits,  na  =  V  (the  velocity  of  the  satellite)  and  since  e  «  1 


(6.57)  becomes 


^e/dt «  '/V  [  sin  v  f,  +  (  2  cos  v  +  e  -  e  co;>-  v)  f2]  (6.58) 


To  determine  the  value  of  fj  and  f2  the  following  assumptions  are  made: 

(1)  the  distance  between  the  Sun  and  satellite  is  infinite;  the  parallax  of  the 
Sun  is  therefore  negligible. 

(2)  the  solar  flux  is  constant  along  the  orbit  when  there  is  no  shadow;  it  is 
time  independent. 

(3)  There  is  no  re-radiation  from  the  surface  of  the  Earth. 


The  pressure,  P,  exerted  by  the  electromagnetic  radiation  from  the  Sun  is 
approximately  4.5  x  10  6  newton/m2  [37].  The  electromagnetic  radiation  pressure 


exerts  a  force  on  the  satellite  proportional  to  its  cross-section/mass  ratio.  Accurate 
modeling  of  the  solar  radiation  acceleration  is  very  difficult  because  the  cross-section 
exposure  to  sunlight  varies  along  the  orbit.  In  addition,  the  acceleration  ceases  each 
time  the  satellite  passes  through  the  shadow  of  the  Earth. 

A  convenient  parameter  to  use  is  the  effective  cross-section/mass  ratio 
defined  by 

o  =  A/m  (1+  e )  typical  values  are  of  the  magnitude  10"2m2/kg 

where 

A  =  cross-sectional  area 
m  =  mass 

€  =  a  factor  determined  from  empirical  data 

Therefore, 

» 

f[  =  -Po  cos  s  sin  i 
f2  =  Pg  sin  s  sin  i 

where  s  =  the  sidereal  angle  of  the  Sun 

i  =  angle  between  the  direction  of  the  Sun 
and  the  orbit  normal 

Equation  (6.58)  becomes 

^e/dt  =  Po/V  [-sin  v  cos  s  sin  i  +  (2  cos  v  +  e  -  e  cos2v)  sin  s  sin  i  ] 

(6.59) 

the  factor  P/V  is  approximately  1 .298  x  lCh4  kg/m2 

To  use  (6.59)  the  true  anomaly  at  which  the  satellite  enters  and  leaves  the 
Earth's  shadow  must  be  known.  The  change  in  eccentricity  over  a  period  of  time  can 
then  be  obtained  by  numerical  integration  of  the  orbit.  For  the  Solar  Power  Satellite 
with  an  area  to  weight  ratio  of  1 .73  m2/kg  the  max  change  in  eccentricity  over  a  year 


varies  between  .025  and  .042  depending  upon  inclination  angle  and  the  ascending 
node  (see  reference  75  for  detailed  figures  and  results).  The  motion  of  eccentricity  is 
periodic  with  a  period  of  one  year. 

To  compensate  for  the  change  in  eccentricity  a  variety  of  bum  cycles  can  be 
used.  If  judiciously  selected  the  bums  can  be  accomplished  in  conjunction  with  the 
tangential  bums  (see  section  6.6)  performed  for  longitude  station  keeping. 
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Lunar  Gravitational  Effect  on  Near-Circular  Orbits 
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